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Final  Report  - Investigation  of 
Finite  Elements  for  Strongly  Nonlinear  Problems 

This  final  report  summarizes  the  results  of  an  evaluation  of  two  new  plate  and 
shell  elements  for  use  in  large  deflection  nonlinear  analysis.  The  elements 
are  of  the  "stability  element"  type,  in  which  nonlinear  strains  are  included  in 
calculations,  with  their  values  optimized  by  added  membrane  displacement  func- 
tions and  special  types  of  elemental  level  constraints.  The  report  summarizes 
the  formulation  and  implementation  of  the  elements,  and  discusses  numerical 
results  in  detail.  Conclusions  are  drawn  regarding  the  effectiveness  of  these 
elements  for  solving  highly  nonlinear  problems  and  also  for  solving  certain 
types  of  linear  shell  structure  problems.  The  crucial  role  of  stepping/ 
iterative  solution  procedures  is  discussed,  and  solution  procedures  of  improved 
convergence  are  described.  Finally  nonlinear  finite  element  alternatives  are 
discussed  from  the  points  of  view  of  accuracy,  computing  cost,  element  size, 
discretization  considerations,  and  solution  convergence. 
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This  report  discusses  the  development  and  evaluation  of  two  new  finite  elements 
for  nonlinear  shell  analysis.  The  elements  are  of  a new  type  called  "stability 
elements,"  specifically  adapted  to  handle  linear  shell  analysis,  and  nonlinear 
effects  in  plates  and  shells  due  to  large  deflections.  The  original  work  on 
these  elements  was  under  NASA-MSFC  Contract  NAS8-30626,  monitored  by  Dr.  John  Key, 
and  accomplished  theoretical  development  and  a partial  computer  code  develop- 
ment. The  present  research  is  under  AFOSR  Contract  F44620-76-C-0130,  monitored 
by  Mr.  William  Walker  and  Lt.  Col.  Joseph  Morgan.  Under  the  present  contract 
the  computer  codes  have  been  completed,  critical  evaluations  made  of  the  perform- 
ance of  the  elements,  and  recommendations  developed  for  further  research. 

The  numerical  performance  of  the  stability  elements  has  been  very  good,  and  it 
appears  clear  that  they  have  a significant  advantage  over  conventional  elements 
for  solving  nonlinear  problems  without  resorting  to  very  small  elements.  The 
present  research  has  defined  modified  and  simplified  formulations  for  the 
elements  which  should  improve  their  performance  still  further.  In  addition, 
the  numerical  work  has  indicated  that  improved  nonlinear  solution  procedures 
should  be  used  in  conjunction  with  the  new  elements,  in  order  to  achieve  the 
large  step  sizes  which  these  elements  can  handle.  Recommendations  are  made  for 
a combination  of  nonlinear,  stability-type  elements  with  nonlinear-step  solution 
procedures. 

Computing  costs  with  the  new  elements  have  been  quite  high.  However,  high 
computing  costs  are  generally  acceptable  in  stepping  solutions  of  highly  non- 
linear problems.  It  appears  that  with  the  use  of  conventional  elements,  in 
smaller  sizes,  to  achieve  comparable  accuracy,  the  costs  would  be  higher  still. 

The  use  of  the  recommended  nonlinear-step  solution  procedure,  in  conjunction 
with  the  stability  elements,  should  decrease  computing  costs  significantly 
through  allowing  larger  step  sizes. 

This  report  describes  the  technical  formulation  of  the  stability  elements  and 
outlines  the  overall  computational  procedures  required  for  their  use.  Numerical 
calculations  are  discussed  for  several  problems  in  order  to  illustrate  the 
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elements'  performance  in  critical  areas.  Solution  procedure  difficulties  are 
discussed,  and  a special,  quasi-nonl inear  solution  procedure  which  was  success- 
ful in  achieving  convergence  with  relatively  large  step  sizes  is  described. 
Emphasis  in  the  report  is  placed  on  the  quadrilateral  element,  since  the  work 
has  indicated  that  it  is  the  better  of  the  two  elements  studied.  There  are 
alternative  approaches  to  the  stability  elements  which  might  be  used  to  solve 
highly  nonlinear  problems  with  varying  degrees  of  success  and  cost.  These  are 
briefly  discussed  in  Section  2.6.  Finally,  conclusions  and  recommendations  for 
further  research  are  discussed  in  Section  3.0. 


2.0  TECHNICAL  DISCUSSION 


This  section  discusses  the  conceptual  basis,  theoretical  development,  computer 
code  arrangement,  compatible  stepping  solution  procedures,  and  numerical  per- 
formance of  the  new  stability  elements.  This  comprises  an  overall  description 
of  the  element  development  from  conception  through  the  evaluations  and  exten- 
sions accomplished  on  the  present  contract.  This  broad  coverage  will  describe 
the  elements  sufficiently  that  reference  to  earlier  documentation  for  detail 
should  not  be  required.  However,  should  details  such  as  equations  or  more  in- 
depth  discussions  be  required,  the  reader  is  referred  to  the  original  contract 
document  (Reference  1)  submitted  to  NASA-MSFC  (NAS8- 30626,  Volume  II). 

The  evaluation  of  the  numerical  performance  of  the  elements.  Section  2.5, 
includes  discussions  of  specific  areas  in  which  their  formulation  requires 
improvement.  These  results  are  the  basis  of  the  recommendations  made  in 
Section  3.0. 

2.1  Description  of  Stability  Elements 

It  has  been  found  that  conventional  finite  elements  may  suffer  serious  loss  of 
accuracy  in  application  to  large  deflection  problems  which  are  strongly  non- 
linear. The  particular  difficulty  of  concern  here  occurs  when  the  finite 
element  formulation  specifically  employs  nonlinear  strain-displacement  equations 
and  thus  may  generate  large  strain  levels  through  nonlinear  behavior.  The 
problem  was  originally  discussed  by  Mallett  in  Reference  2,  by  Haftka,  Mallett, 
and  Nachbar  in  Reference  3,  and  by  Berke  and  Mallett  in  Reference  4.  The  work 
in  Reference  3 suggested  a procedure  for  resolving  the  difficulty  for  beam 
elements.  The  present  work  has  extended  this  procedure  to  the  two  dimensional 
case  of  plates  and  doubly  curved  shells,  and  refers  to  the  elements  as  HMN 
elements,  in  recognition  of  the  authors  of  Reference  3. 

The  basic  difficulty  occurs  because  the  nonlinear  strain-displacement  equations 
contain  squares  and  products  of  the  displacement  gradients,  in  particular,  the 
bending  slopes  of  the  midsurface  of  the  plate  or  shell.  Since  the  bending 
displacements  of  such  elements  conventionally  are  second,  third,  or  higher 
degree  polynomials,  the  squares  and  products  of  their  derivatives  are  of  degree 


two,  three,  four,  or  higher.  The  membrane  strains  therefore  contain  terms,  due 
to  nonlinear  behavior,  which  are  of  a higher  polynomial  degree  than  those  which 
result  from  the  derivatives  of  the  membrane  displacements  themselves.  These 
higher  degree  polynomial  forms  also  occur  in  shell  elements  in  linear  analysis, 
due  to  the  presence  of  initial  curvature.  The  high  degree  membrane  strains 
induced  by  nonlinearity  or  initial  curvature  cannot  be  •erased"  by  the  lower 
degree  element  membrane  displacement  function  derivatives.  The  result  is 
excessive  strain  energy,  and  hence  stiffness,  of  the  finite  element  representa- 
tion. In  contrast,  in  actual  physical  behavior,  plate  and  shell  structures 
behave  in  such  a way  as  to  minimize  the  participation  of  higher  degree  deforma- 
tion forms,  since  such  forms  would  excessively  absorb  strain  energy  and  thus 
would  cause  excessively  high  stiffness  in  structural  behavior.  The  actual 
physical  action  takes  place  through  small  in-surface  adjustments  of  the  membrane 
displacements,  which  naturally  work  to  reduce  the  strain  energy  by  eliminating 
unnecessarily  complex  strain  states.  The  membrane  displacements  which  accom- 
plish this  are  necessarily  of  a rapidly  varying,  i.e.,  high  polynomial  degree, 
type. 

Stability  Element  Formulation 

The  approach  of  Reference  3 was  to  introduce  axial  displacement  forms  for  the 
beam  element  through  the  5th  degree  polynomial  forms.  This  provided  a quartic 
axial  strain  which  was  used  to  ■erase"  the  quartic  nonlinear  axial  strain 
resulting  from  the  square  of  the  bending  slope.  The  method  employed  to  deter- 
mine the  amplitudes  of  the  added  axial  displacement  forms  was  minimization  of 
the  potential  energy  on  the  elemental  level.  This  corresponds  to  the  physical 
behavior  by  which  the  structure  seeks  a minimum  energy  state.  The  authors  also 
showed  that  the  same  result  can  be  obtained  by  directly  constraining  the  axial 
strain  to  eliminate  its  high  degree  polynomial  components. 

The  procedures  of  Reference  3 appear  to  be  extendable  to  apply  to  plate  and 
shell  elements.  The  basic  approach  is  to  start  with  an  available  element  which 
has  been  used  for  linear  analysis,  to  add  supplementary  membrane  (HMN)  functions 
to  control  the  high  degree  nonlinear  membrane  strains,  and  to  derive  appropriate 
constraint  equations  to  effect  this  control.  This  was  don§  with  the  BCIZ  and 


AZI  (References  5,  6)  elements.  Considerable  additional  complexity  occurs,  for 
several  reasons:  there  are  three  membrane  strains  to  control  and  two  displace- 
ment components  to  deal  with,  rather  than  the  single  axial  strain  and  displace- 
ment of  the  beam  problem;  constraints  on  the  added  displacement  forms  must 
avoid  creating  inter-element  incompatibilities;  the  strain  constraints  must  be 
formulated  in  two  dimensions,  and  the  equations  are  difficult  to  deduce.  The 
overall  procedures  required  to  deal  with  the  plate  and  shell  elements  is 
described  in  detail  in  Reference  1.  Briefly,  it  is  as  described  below.  The 
development  is  based  on  extensions  of  the  triangular  BCIZ  (Reference  5)  and 
quadrilateral  AZI  (Reference  6)  elements.  Double  curvature  and  fully  coupled 
membrane  and  bending  displacement  states  are  included. 

o Higher  degree  polynomial  forms  in  the  membrane  strain  equations,  result- 
ing from  large  deflection-induced  nonlinearities,  and  also  from  initial 
curvature,  for  shell  elements,  are  determined.  Thus,  amplitudes,  in  terms 
of  nodal  bending  freedoms,  of  polynomial  forms  such  as  x2y,  xy2,  x3,  y3, 
x2y2,  etc.,  are  determined. 

o Supplementary  membrane  displacement  forms  are  deduced  such  that  the  added 
forms  are  able  to  produce  these  same  polynomial  strain  terms.  The  supple- 
mentary membrane  functions  must  form  a complete  set,  in  combination  with 
the  basic  membrane  functions  of  the  element.  The  high  degree  polynomial 
strain  terms  must  be  produced  as  independent  functions. 

o A set  of  strain  constraints  is  developed  to  reduce  the  polynomial  degree 
of  the  membrane  strains  to  the  same  degree  which  they  have  in  the  basic 
element  formulation.  The  constraints  must  avoid  inter-element  incompati- 
bility while  providing  as  complete  as  possible  a control  over  the  high 
degree  membrane  strains.  Dependencies  among  the  strain  constraint  condi- 
tions must  be  avoided;  this  can  be  a troublesome  problem. 

o Those  supplementary  (HMN)  membrane  functions  which  do  not  participate  in 
the  above  constraints,  but  which  are  present  in  order  to  keep  a complete 
set  of  functions,  are  reduced  out  by  potential  energy  minimization  at  the 
elemental  level.  The  minimum  energy  constraint  cannot  be  used  for  HMN 
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functions  at  inter-element  boundaries,  because  inter-element  displacement 
incompatibilities  would  result. 

In  applications,  the  above  has  led  to  procedures  in  which  the  added  displace- 
ments include  high  degree  forms  defined  independently  along  the  sides  of  the 
elements  and  in  their  interiors,  and  in  which  the  constraints  include  both 
explicit  strain  constraints  and  minimum  potential  energy  constraints,  all  on 
the  elemental  level.  When  the  added  displacements  are  parallel  to  a particular 
line  (say,  a side)  in  the  element,  the  amplitudes  of  the  added  displacements 
are  completely  fixed  by  the  values  of  the  bending  displacement  along  that  line, 
with  numerical  values  defined  by  the  constraint  equations.  This  is  the  case 
for  most  of  the  HMN  constraints  used  in  the  two  elements,  and  assures  inter- 
element compatibility  of  displacements  parallel  to  the  sides  of  the  elements. 
For  the  triangular  element,  HMN  functions  normal  to  the  sides  of  the  elements 
were  found  to  cause  inter-element  incompatibilities,  and  hence  were  dropped 
from  consideration.  In  the  numerical  work  to  date  with  the  quadrilateral 
element,  supplementary  membrane  displacements  normal  to  the  sides  of  the  ele- 
ments have  in  some  cases  been  used.  These  functions  serve  to  eliminate  high 
degree  nonlinear  membrane  shear  strains.  The  element  derivation  contains 
such  normal -to-the-side  displacements  of  cubic  and  quartic  variations.  It  has 
been  reasoned  (Reference  1)  that  the  cubic  form  would  produce  unacceptable 
inter-element  incompatibility,  and  these  functions  are  not  retained.  The 
quartic  form  was  felt  to  be  acceptable,  and  is  currently  retained  in  calcula- 
tions. The  inter-element  incompatibility  associated  with  the  quartic  function 
is  nonzero  but  generally  very  small. 

Figure  1 illustrates  the  manner  in  which  the  explicit  strain  constraints  are 
imposed  for  the  two  elements  in  the  computer  programs.  For  the  triangle,  the 
extensional  strain  parallel  to  the  three  sides,  designated  by  e$,  is  controlled 
as  a polynomial  function  of  the  side-length  coordinate,  s.  This  control  is 
through  the  4th  degree  polynomial.  For  the  quadrilateral,  the  extensional 
strain  parallel  to  the  edges  and  the  mid-line  of  the  element,  in  both  the  x and 
y directions,  are  controlled  through  the  2nd  degree  polynomial  form.  Thus,  ex 
is  controlled  along  lines  1-3,  8-4,  7-5,  and  ey  along  lines  1-7,  2-6,  and  3-5 
in  the  figure,  using,  respectively,  cubic  functions  for  u and  v.  In  addition, 
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the  shear  strain  exy  is  controlled  for  the  cubic  polynomial  form  along  the  four 
sides  of  the  element.  This  is  accomplished  through  added  displacements  of 
quartic  form,  normal  to  the  sides  of  the  elements,  as  discussed  above. 

All  HMN  constraints  are  formulated  for  the  parent  element  of  the  isoparametric 
family.  . This  is  essential  because  the  displacement  functions  are  expressed  in 
terms  of  the  parent  element  coordinates.  It  is  essential  to  account  for  the 
general  isoparametric  element  in  the  HMN  constraints.  This  is  not  conveniently 
done  through  the  normal  Jacobian  transformation  of  the  isoparametric  family, 
since  distinct  polynomial  forms,  rather  than  numerical  values  at  integration 
points,  are  required  for  the  HMN  constraints.  Consequently,  it  was  necessary 
to  re-derive  the  HMN  constraints,  dealing  at  the  outset  with  the  general  iso- 
parametric element.  This  was  done  through  the  use  of  element  side  and  mid-line 
(see  Figure  2)  length  scale  factors  and  similar  scale  factors  relating  to 
differentiations. 


I 

I 

I 

I 


Coordinate  Systems  and  Updating 

The  application  of  the  HMN  strain  constraint  procedure  and  the  intended  purposes 
of  the  new  elements  have  implications  regarding  coordinate  systems  and  stepwise 
updating.  The  elements  are  to  handle  large  rotation  problems  and  are  to  avoid 
cumulative  error  by  computing  total  nonlinear  strains  rather  than  summing 
increments.  These  requirements  are  best  satisfied  by  using  element  coordinate 
systems  which  follow  the  elements  throughout  the  deformation.  If  a fixed 
system,  such  as  a fixed  cartesian  or  fixed  shell  surface  intrinsic  system  were 
used,  and  the  stepwise  calculation  implemented  by  updating  the  element  position 
within  this  system,  a problem  of  role  exchange  between  the  element  displacement 
forms  would  occur.  In  this  role  exchange,  the  displacement  forms  (e.g.,  cubics) 
intended  for  the  bending  displacements,  and,  say,  directed  along  a locally 
fixed  z-axis,  would  after  a large  rotation  be  partly  directed  along  material 
lines  lying  in  the  plane  of  the  deformed  element.  Conversely,  the  membrane 
displacements  along  locally  fixed  x and  y axes  would,  after  a large  rotation, 
be  partly  directed  along  a material  line  normal  to  the  deformed  element.  This 
difficulty  limits  the  allowable  magnitudes  of  the  rotations  for  elements  with 
distinctly  different  membrane  and  bending  displacement  forms,  in  formulations 
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using  fixed  coordinate  systems  for  displacement  definitions.  It  suggests  that 
"following"  coordinate  systems  with  updating  be  used  for  nonlinear  elements. 

A second  factor  reinforces  this  conclusion  in  the  case  of  the  stability  ele- 
ment approach.  The  nonlinear  strain  equation  for  ex,  referred  to  locally 
cartesian  coordinates,  is 

‘x  - “-X  ♦ 1/2 

2 2 2 
The  term  w,x  is  important  even  for  small  rotations;  u,x  and  v>x  become 

important  terms  for  large  rotations.  As  discussed  earlier,  the  membrane 
functions  u and  v contain  high  degree  polynomial  forms,  to  implement  the 
HMN  methodology.  If  the  underlined  terms  in  the  above  equations  are  retained, 
these  high  degree  components  of  u and  v will  create  undesirably  high  degree 
polynomials  in  ex.  This  is  precisely  the  effect  which  the  HMN  procedure  is 
intended  to  avoid.  Consequently,  in  formulating  the  HMN  elements,  it  appears 
to  be  necessary  to  limit  the  magnitudes  of  the  rotations  through  the  use  of 
coordinate  systems  which  follow  the  elements  throughout  the  entire  deformation 
process,  and  correspondingly  to  omit  terms  in  the  strain  displacement  equa- 
tions of  the  type  underlined  above.  A system  of  cartesian,  updated,  element 
baseplane  coordinate  systems  was  developed  for  this  purpose.  The  rotations 
which  occur  during  a given  load  step  are  small,  and  after  the  step  has  been 
computed,  and  the  cartesian  baseplane  coordinate  systems  updated,  the  total 
rotations  of  the  material  elements  relative  to  their  updated  coordinate  sys- 
tems remain  small.  Thus,  total  strains  can  be  computed  accurately  with  all 
nonlinear  terms  of  the  types  underlined  above  omitted  from  the  membrane  strain 
equations.  In  use  for  finite  element  analysis  with  reasonably  small  elements, 
this  approach  guarantees  that  element  material  slopes  will  always  be  small 
relative  to  the  baseplane  of  reference.  Hence,  it  is  possible  to  use  a shal- 
low shell  type  of  deformation  definition,  with  its  attendant  simplification  of 
formulation.  This  is  a significant  gain  for  the  nonlinear  elements  under  con- 
sideration, since  non  linearities  develop  through  accumulated  element  slopes 
rather  than  through  cumulative  changes  of  the  shell  geometrical  parameters, 
referred  to,  say,  a lines-of-curvature  coordinate  system. 
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The  baseplane  following  coordinate  system  procedure  is  illustrated  by  Figure  2. 
In  this  figure,  coordinate  axes  (xQ,  zq)  are  used  to  set  up  the  equations  for 
the  first  load  step,  which  terminates  with  the  element  in  the  position 
described  by  the  updated  axis  system  (x^,  z-j).  Similarly,  load  steps  two  and 
three  cause  displacements  to  the  positions  described  by  the  updated  systems 
(x2»  *2^  ancl  (x3»  z3)»  etc*  *n  updating  the  coordinate  systems  in  this  way, 
the  total  displacements  of  the  element  are  transformed  to  the  new  axes,  with 
only  the  rigid  motions  removed.  Used  in  conjunction  with  a Lagrangian  strain 
definition,  this  approach  is  best  termed  a "start-over-total -Lagrangian" 
approach,  since  total  strains  are  repeatedly  computed  from  new  starting  orien- 
tations. The  displacement  transformation  is  not  a straightforward  geometrical 
transformation  of  displacement  increments  in  the  usual  sense.  Instead,  it  is 
a special  transformation  based  on  the  idea  of  computing  total  element  displace- 
ments, referenced  to  the  updated  coordinate  system,  measured  between  the 
deformed  element  position  and  an  undeformed  element  suitably  mapped  onto  the 
updated  system.  This  procedure  assures  that  total  strains  can  be  correctly 
computed  using  the  transformed  displacements  referred  to  any  updated  coordi- 
nate system.  Figure  3 illustrates  the  basis  of  the  transformation  for  the 
simple  case  of  a beam  element.  An  element  is  shown  referred  to  an  initial 
coordinate  system  (xq,  zq)  and,  after  a deformation  step,  to  an  updated  system 
(x-j,  z.|).  The  undeformed  shape  of  the  element  is  shown  on  both  coordinate 
systems;  the  element  end  points  always  lie  on  the  x-axis.  Displacements  are 
defined  to  be  the  vector  differences  between  points  on  the  deformed  and 
undeformed  elements,  referred  to  the  particular  coordinate  system  in  question. 
This  definition  suffices  for  the  calculation  of  strain.  The  element  in  system 
(xQ,  zq)  has  accumulated  displacements  which,  through  prior  transformations, 
are  referred  to  this  system.  Thus,  for  a point  P,  the  accumulated  displace- 
ments can  be  denoted  by  the  vector  notation 


(V  "op> 

In  which  the  subscripts  denote  both  the  reference  system  and  the  point  P.  The 
incremental  displacements  which  occur  during  the  present  step  are  also  referred 
to  the  (xQ,  zQ)  system,  and  are 
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The  position  vector  of  point  P on  the  undeformed  element  in  the  (xQ,  zQ)  sys- 
tem is  denoted  by 


<V  V 


Thus  the  point  P after  the  current  step  has  the  position 


(X0P-  V + (V  “op1  * (4V  “op1 


Referred  to  the  (x^,  z^)  system,  the  position  of  point  P on  the  undeformed 
element  is  denoted  by 


(v  v 

It  is  noted  that  the  individual  components  satisfy  x,  = x and  z,  = z , 

lp  op  lp  op 

since  these  values  refer  to  the  undeformed  state.  The  total  accumulated  dis- 
placements at  the  end  of  the  step,  referred  to  the  (x^,  z^)  system  are 

*ulp’  wlp} 

The  definitions  of  these  components  are  as  shown  on  Figure  3.  The  position  of 
point  P can  now  be  given  in  both  coordinate  systems.  Equating  these  necessarily 
identical  vector  quantities  gives 


(x„  . z ) + (u , w)  + (au , Aw  ) 

op  op  op  op  op  op 


R * <V  V * (V  "lp* 


where  R is  the  translation  of  the  origin.  This  equation  is  solved  for  (ulp, 
w-|p),  in  component  form,  by  employing  the  appropriate  unit  vectors  of  the  two 
coordinate  systems  and  forming  dot-products  with  both  sides  of  the  equation. 
Differentation  of  the  vector  components  ulp  and  w^p,  with  respect  to  the 
updated  baseplane  coordinates,  treated  as  material  coordinates,  yields  trans- 
formed displacement  derivative  freedoms,  if  required,  and  also  all  of  the 
quantities  required  for  calculation  of  the  total  nonlinear  strains.  A detailed 


set  of  equations  of  this  transformation  procedure  is  given  in  Reference  1.  It 
is  seen  that  this  system  of  treating  displacements  handles  the  element  displace- 
ment functions  as  displacements  in  the  directions  of  the  convected  baseplane 
coordinates.  Each  incremental  displacement  consists  of  element  displacement 
function  increments  referred  to  an  updated  element  baseplane. 

There  is  an  equivalent  alternative  to  the  above  approach  for  the  direct  calcu- 
lation of  purely  nodal  displacements.  This  is  to  maintain  cumulative  nodal 
displacement  values  referred  to  a fixed  global  coordinate  system,  and  use 
these  values  to  transform  back  to  the  successively  updated  element  systems. 

Such  a procedure  is  a useful  part  of  the  solution  procedure  in  any  case, 
because  the  global  values  are  a convenient  means  of  determining  the  orienta- 
tions of  the  successive  element  baseplanes.  However,  this  procedure  does  not 
provide  nodal  displacement  derivative  freedoms,  which  are  needed  for  some  types 
of  elements.  Both  methods  are  used,  as  appropriate,  in  the  computer  coding  of 
the  two  elements  under  study  in  this  contract.  The  triangular  element  (based 
on  the  BCIZ  element,  Reference  5)  uses  nodal  derivative  freedoms.  For  this 
element,  the  above  described  transformation  was  found  essential  for  transform- 
ing the  derivative  freedoms.  The  quadrilateral  element  (based  on  the  AZI 
element.  Reference  6)  uses  nodal  displacements  and  nodal  rotations  as  free- 
doms. For  this  element,  since  rotations  rather  than  derivative  freedoms  are 
used,  a different  type  of  transformation  was  required.  This  particular  trans- 
formation was  developed  during  the  checkout  of  the  quadrilateral  element,  and 
is  outlined  below. 

Figure  4 shows  a portion  of  a deformed  element  of  the  AZI  type,  in  which  trans- 
verse shear  deformations  are  permitted.  For  simplicity  a beam  element  is  con- 
sidered here.  The  displacement  state  of  the  element  is  referred  to  its  succes- 
sive baseplanes.  The  start-of-step  baseplane  is  denoted  here  by  the  xQ  axis, 
and  the  end-of-step  baseplane  by  the  x-j  axis.  The  displacement  quantities  for 
the  start-  and  end-of-step  states,  for  a point  P on  the  midline  of  the  element 
and  for  the  fiber  PS,  defined  to  be  the  original  normal  to  the  midline  in  the 
undeformed  state,  are  indicated  in  the  figure.  The  subscripts  refer  to  the 
coordinate  system  to  which  the  displacements  are  referred.  The  transformation 
described  earlier  in  this  section  permits  the  calculation  of  w^  from  the  values 
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w + Aw  . Here  we  are  concerned  with  the  determination  of  e,  from  the 
op  op  lp 

values  of  eQp  and  A9Qp.  A procedure  which  immediately  suggests  itself  to 
make  this  determination  is  a simple  subtraction  of  angles,  including  the  angle 
change  between  the  successive  baseplanes.  This  is  not  the  procedure  used, 
however.  Instead,  the  method  developed  is  based  on  representing  the  material 
vector  PS  in  the  successive  coordinate  systems  and  deducing  the  needed  angles 
from  its  components.  This  method  is  convenient  in  the  three  dimensional  case, 
and  also  incorporates  the  effects  of  the  rotations  on  the  z-variation  of  the 
bending  displacement,  w.  This  was  found  essential  for  accurate  calculation 
of  nonlinear  strains.  The  transform  step  is  briefly  described  here,  since  it 
is  not  considered  in  Reference  1.  The  method  described  here  is  used  in  the 
compute*'  code  for  the  quadrilateral  element.  A more  exact  procedure,  which 
would  be  required  for  larger  single  step  displacements  or  for  shell  elements 
joining  non-tangential ly , is  described  in  Section  2.4.  We  assume  as  before 
that  all  local  and  single  step  angles  and  strains  are  small  compared  to  unity 
and  that  e and  sin  e can  be  taken  as  equal.  Then  the  vector  PS  in  the  xQ 
system  is  given  by 


where  | PS | is  the  original  length  of  PS,  subscripts  x and  y denote  rotations 

about  the  x and  y axes,  and  the  vectors  "\T  are  unit  vectors  in  the  x^  system, 
o o o 

The  assumption  that  the  incremental  angles  of  rotation  are  sufficiently  small 
to  be  added  in  any  order  is  implicit  in  this  equation.  The  third  component 
of  the  row  vector  retains  the  squared  terms  in  order  to  included  the  effect 
of  foreshortening  of  the  normal  due  to  its  rotations.  For  simplicity  we 
abbreviate  this  equation  as  follows 

PS  = |PS|  • L e0  J • {vQ\ 

After  the  deformations  of  the  current  step, 

PS  + APS  = | PS | • L eQ  + A60  J • {^T| 
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| PS | is  unchanged  because  the  condition  of  inextensible  normals  is  imposed. 
This  condition,  and  its  use  in  determining  the  nonlinear  transverse  shear 
strain,  discussed  later,  are  the  reasons  for  retaining  the  two  squared  terms 

. Vl-'e*' 


in  the  zQ  component 


-e~  -e^  ' 
opy  opx 


The  vector  PS  + APS  is  easily  transformed  to  the  x^  system  by  means  of  a 

direction  cosine  matrix  of  the  cartesian  coordinate  transformation  x — ► x, 

o 1 

Denoting  this  matrix  by  [ X ] , there  results 

PS  + APS  = L 0Q  + A0q  J * [x]  = L e]  J • { v*j  | 

The  definition  of  [xj  consistent  with  the  above  equation  is 


Lx]  - 


in  which  the  indicated  operations  are  scalar  products  and  the  subscripts  of 
the  unit  vectors  refer,  respectively,  to  the  coordinate  system  to  which  they 
belong,  and  their  particular  component  direction,  in  the  order  x,  y,  z.  The 
updated  vectorjcj,  is  thus  given  by 
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and  the  needed  values  0,  and  e,  are  obtained  as  the  first  component  and 

i py  i px 

the  negative  of  the  second  component  of  {e^}.  The  matrix  [x]  is  readily 
obtained  in  the  computer  code  from  the  product  of  the  transformations  between 
the  element  baseplane  coordinate  systems  and  the  global  coordinate  system, 
for  the  start-of-step  and  end-of-step  conditions. 

The  bending  displacements,  w,  are  also  transformed  between  the  start-of-step 
and  end-of-step  systems.  The  result  of  this  transformation,  in  conjunction 
with  the  known  displacement  shape  functions  of  the  element,  permits  calcula- 
tion of  the  bending  slope  quantities  w,  and  w,  , referred  to  the  end-of-step 

x y 
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system.  In  doing  this,  the  original  element  dimensions  and  shape  are  used, 

consistent  with  the  approach  of  using  material  coordinates  and  obtaining  a 

nonlinear  strain  of  the  Lagrangian  type.  The  element  w-displacement  functions 

are  considered  to  be  convected  to  the  updated  baseplane  coordinate  system,  in 

order  to  determine  the  derivatives  w,  and  w,  . The  end-of-step  values  of  w,  , 

a y x 

w,  , and  the  rotations,  0 and  6 , suffice  to  obtain  the  total  nonlinear 

J \ja  py 

transverse  shear  strain.  Further  details  of  this  calculation  are  discussed 
below. 


Shell  Strain-Displacement  Equations 


A final  item  will  complete  this  description  of  the  conceptual  basis  of  the  two 
shell  elements  under  study.  It  was  noted  earlier  in  this  section  that,  to 
obtain  a workable  theory  for  stepping  out  nonlinear  solutions,  it  is  preferable 
to  update  the  element  baseplane  in  order  to  maintain  a small  angle  relationship 
between  the  element  and  its  reference  system.  The  resulting  availability  of 
the  updated  baseplane  strongly  suggests  the  use  of  a shallow  shell  theory  in 
which  the  displacements  and  forces  are  referred  to  the  cartesian  baseplane 
system  rather  than  to  the  shell  intrinsic  curvilinear  system.  The  simplifica- 
tions obtained  in  this  way  are  particularly  helpful  for  nonlinear  problems. 
Consequently,  it  was  decided  to  base  the  element  development  on  a shallow  shell 
type  of  theory,  using  cartesian  displacement  definitions,  with  baseplane 
updating  to  permit  extension  to  the  large  deflection  regime.  Since  shallow 
shell  theories  do  not  as  a rule  use  purely  cartesian  displacement  definitions, 
it  was  necessary  to  derive  a new  set  of  equations,  starting  from  basic  princi- 
ples. To  do  this,  a tensor  approach  was  used,  and  the  strains  were  derived 
from  the  changes  of  the  metric  tensor  between  the  undeformed  and  the  deformed 
states.  This  is  a superior  approach,  which  naturally  provides  all  of  the 
nonlinear  effects  of  large  displacements,  without  recourse  to  the  more  conven- 
tional geometric-deductive  procedures.  The  derivations  for  the  triangular 
(Kirchhoff-type)  and  quadrilational  (transverse  shear  strain  included)  elements 
are  given  in  Reference  1.  Numerical  work  has  pointed  to  a very  important 
result  of  this  overall  procedure.  The  use  of  cartesian  displacements,  specifi- 
cally the  w-displacement,  leads  to  a much  different  definition  of  the  direct 
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! membrane  strains  than  the  conventional  theory  (both  deep  and  shallow  shells)  in 

which  w is  normal  to  the  shell  midsurface.  Conventionally,  the  membrane  strain 
includes  terms  such  as  w/R,  where  R is  a shell  midsurface  radius  and  w is 
normal  to  the  midsurface.  This  term  is  replaced,  when  cartesian  displacements 
are  used,  by  terms  like  w,°  • w,x,  where  w°  is  the  initial  w-shape  of  the  shell 
element,  referred  to  its  baseplane,  and  the  commas  denote  differentiation.  For 
a given  w distribution  over  the  element,  these  two  terms  provide  very  different 
forms  of  contributions  to  the  membrane  strain.  They  are  made  equivalent  by 
differing  forms  of  the  membrane  displacement,  u,  which  contributes  the  term  u,x 
to  the  strain.  It  can  be  reasoned  that,  for  conventional  strain  equations, 
used  in  a finite  element  application,  the  membrane  displacements  are  called 
upon  in  part  to  supplement  the  bending  displacement  in  order  to  achieve  satis- 
factorily strain-free  rigid  body  motions  of  the  elements.  This  generally 
requires  very  competent  membrane  displacements,  particularly  for  low  energy 
deformations  such  as  inextensional  bending.  For  the  cartesian-based  strain 
equations,  rigid  motions  are  automatically  obtained,  and  the  membrane  displace- 
ments are  called  upon  to  supplement  the  bending  displacement  in  order  to  achieve 
smoothly  varying  membrane  strain  states.  In  both  cases,  the  use  of  membrane 
displacements  of  higher  polynomial  degree  than  the  bending  displacements  is 
required  for  accurate  solution  of  a wide  class  of  shell  problems.  All  of  the 
above  pertains  to  the  case  of  linear  analysis,  and  is  governed  primarily  by  the 
effects  of  shell  curvature.  For  nonlinear  analysis,  a similar  situation  exists 
except  that  the  membrane  functions  must  be  sufficiently  competent  to  handle 
large  rigid  motions  and  changes  in  shell  geometry  (curvatures,  twist)  in  their 
efforts  to  achieve  stress-free  rigid  motions  and  smoothly  varying  strain  states. 
It  is  noted  at  this  point  that  the  HMN  functions,  which  increase  element  mem- 
brane displacement  competence,  are  beneficial  for  linear  and  nonlinear  analysis 
with  shell  theories  of  either  the  conventional  or  the  cartesian-based  types. 
Example  problem  #1,  Section  2.4,  is  a case  in  which  the  HMN  functions  are 
crucial  to  an  accurate  problem  solution  in  the  linear  case. 

In  approaching  the  derivation  of  the  shallow  shell  type  of  nonlinear  strain 
equations,  consideration  was  given  to  the  results  of  Reference  7,  in  which  it 
was  shown  that  shallow  and  deep  shell  theories  can  give  very  different  results 
for  certain  types  of  problems.  It  was  decided  that  it  is  the  absence  of 
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certain  types  of  forces  in  the  membrane,  i.e.,  baseplane,  equilibrium  equations 
of  conventional  shallow  shell  theory  which  leads  to  significant  errors  in 
particular  problem  solutions  (Reference  7).  The  terms  in  question  are  the 
products  of  the  transverse  shear  forces  with  the  shell  midsurface  slope  angles. 
These  products  provide  membrane-direction  forces,  resulting  from  the  transverse 
shear  forces.  These  are  omitted  in  conventional  shallow  shell  theory  as  a part 
of  the  basic  assumption  that  the  transverse  shear  forces  are  small.  They  can 
be  either  retained  or  omitted  in  the  cartesian-based  theory,  depending  on 
retention  of  certain  small  terms  in  the  strain-displacement  equations.  The 
types  of  problems  in  which  the  error  due  to  this  omission  can  be  large  are 
those  in  which  there  is  strong  bending  combined  with  small  membrane  stress 
levels,  either  over  the  whole  shell  or  in  critical  local  regions.  It  would 
appear  that  nonlinear  large  deflection  behavior  would  accentuate  these  errors, 
for  two  reasons:  large  deflections  are  most  likely  to  occur  in  strong  bending 
rather  than  strong  membrane  deformation  problems;  the  rotations  of  the  large 
deflection  state  will  increase  the  importance  of  those  force  components  which 
are  conventionally  ignored  in  the  shallow  shell  theory.  Consequently,  shell 
equations  were  derived  which  retain  the  simplifications  of  shallow  shell  theory 
without  omitting  these  particular  terms  in  the  equilibrium  equations.  The 
terms  which  were  retained  in  the  strain-displacement  equations  are  in  the 
definitions  of  bending  and  twisting  deformations,  and  involve  products  of 
membrane  displacement  gradients  and  shell  curvatures  and  twist.  Both  initial 
and  subsequent  curvature  and  twist  are  involved,  so  that  the  added  terms  serve 
to  incorporate  nonlinearities  into  the  bending  and  twisting  moments.  Details 
of  this  derivation  and  the  strain-displacement  equations  are  given  in  Reference 
1. 

During  the  numerical  evaluations  of  the  quadrilateral  element,  it  was  found 
that  physically  unexplainable  large  residual  lateral  forces  and  transverse 
shear  strains  occurred  in  nonlinear  bending  problems.  The  cause  of  this 
behavior  was  traced  to  certain  omitted  nonlinear  terms  in  the  transverse  shear 
strains,  as  derived  in  Reference  1.  The  original  derivation  included  the 
products  of  the  rotations  with  the  membrane  displacement  gradients,  but  omitted 
terms  involving  products  of  the  bending  slopes  with  the  z-direction  derivatives 
of  the  bending  displacements,  i.e.,  the  terms  w,  *w,  and  w,  *w,  . These  terms 
were  omitted  on  the  basis  that  all  displacement  z-derivatives  were  dropped  in 
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the  original  derivation,  through  direct  use  of  the  rotations  of  the  normals, 

0 and  0 . However,  it  can  be  easily  shown  that  the  retained  nonlinear  terms, 
x y 

u,*0  ,u,-0  ,v,*0  , and  v,*e  , are  in  most  problems  very  nearly  cancelled  by  the 
xx  yx  xy  yy 

omitted  terms,  w,  • w,  and  w,  • w,  , in  the  nonlinear  strain  equations  for 

fc  A t J 

the  transverse  shear  strains.  Thus,  the  omission  of  the  latter  terms  resulted 
in  large  nonlinear  contributions  to  the  transverse  shear  strains,  producing 
corresponding  large  residual  shear  forces.  Coupling  with  the  bending  slopes 
caused  serious  loss  of  accuracy  and  convergence  difficulties  in  nonlinear 
problems.  When  the  omitted  terms  were  included,  the  transverse  shear  strain 
was  reduced  to  reasonable  values,  and  the  difficulties  with  the  residual  loads 
were  eliminated. 

To  retain  the  effect  of  w,z  in  an  element  which  retains  only  the  freedoms  u,  v, 

w,  0 , and  0 is  a somewhat  difficult  task.  The  value  of  w,  was  constructed 
x y z 

from  the  condition  of  inextensional  normals  and  the  values  of  0 and  e , and 

x y 

has  the  form,  for  small  incremental  rotations, 

w>2  ' - ? (ex  + ey> 

Hence,  the  terms  w,  • w,  and  w,  • w are  cubic  in  the  displacement  magnitudes. 

j L. 

The  basic  strain  formulation  of  the  elements  (Reference  1)  is  of  the  second 
degree  in  the  displacements,  and  does  not  readily  permit  including  the  cubic 
terms.  Consequently,  the  added  nonlinear  terms  were  included  in  the  calcula- 
tion of  strains,  stresses,  and  residual  loads,  but  not  in  the  formation  of  the 
element  stiffness  matrices.  This  type  of  approximation  may  slow  the  conver- 
gence of  stepping/iterative  calculations,  but  does  not  affect  the  accuracy  of 
converged  solutions.  In  future  work  it  would  be  desirable  to  include  the  extra 
terms  at  the  stiffness  matrix  generation  stage  of  the  calculations. 

2.2  Computational  Procedures 

This  section  outlines  the  sequence  of  the  computational  procedures  used  for  the 
two  shell  elements,  with  particular  emphasis  on  those  calculations  which  are 
unique  to  the  stability  elements  and  to  the  coordinate  system  updating  used. 
Details  such  as  equations  and  matrix  definitions  are  given  in  Reference  1. 
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The  purpose  here  is  to  describe  the  overall  scheme  and  magnitude  of  the  compu- 
tational effort. 

Coordinate  Systems 

Each  element  uses  three  principal  coordinate  systems:  the  baseplane  system; 
the  "solution"  systems,  which  are  cartesian  nodal  systems;  and  the  global 
system.  These  are  illustrated  by  Figure  5.  The  solution  systems  are  nodal 
triads  which  are  averages  of  the  joining  element  baseplane  systems.  During 
a single  solution  (load)  step,  a single  solution  system  is  used  without 
updating,  for  each  node  point  of  the  structure.  The  solution  systems  are 
updated  at  the  outset  of  each  solution  step.  Iterative  corrections  to  the 
equilibrium  state  are  made  within  each  solution  step,  based  on  residual  loads 
evaluated  at  each  iteration.  Each  residual  load  evaluation  utilizes  an  updated 
element  baseplane  for  the  evaluation  of  the  strains,  stresses,  and  residuals. 
That  is,  the  displacements  of  the  previous  iteration  are  used  to  update  the 
element  baseplane,  the  total  deformation  state  is  transformed  to  the  updated 
baseplane,  and  total  nonlinear  strains  and  stresses  are  thereby  computed. 

The  virtual  work  integration,  for  the  residual  load  evaluation,  uses  virtual 
displacement  increments  which  are  likewise  referred  to  the  updated  baseplane, 
and  are  increments  from  the  total  deformation  state  referred  to  that  baseplane. 
The  repeated  updating  of  the  baseplane  systems  at  each  iteration  is  done  to 
assure  an  accurate  calculation  of  the  total  Lagrangian  strain,  even  though  the 
displacements  and  rotations  of  the  step  may  be  large.  This  is  probably  not 
necessary  for  well -posed  stepwise  loadings,  and  may  be  a candidate  for  code 
simplification  in  the  future.  The  final  updated  baseplane  systems  of  a given 
load  step,  i.e.,  those  corresponding  to  the  converged  solution  for  the  step, 
are  the  start-of-step  baseplane  systems  for  the  next  load  step.  Figure  6 
illustrates  schematically  the  use  of  these  coordinate  systems  in  a two  step 
problem.  The  figure  also  gives  the  names  of  the  transformation  matrices 
between  the  coordinate  systems,  for  later  reference,  and  indicates  thereby 
which  transformations  require  updating  for  new  load  steps  and  for  iterations 
within  a single  load  step. 
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! For  the  triangular  element  only,  one  other  transformation  of  a coordinate  sys- 

tem type  is  required.  This  transformation  makes  the  transition  between  the 
so-called  "deformational"  freedoms  of  the  original  derivation  of  this  element 
(Reference  5)  and  the  baseplane-referenced  freedoms  which  are  used  in  the 
transformations  [TD]  (see  Figure  6).  It  is  called  "TSTAR.” 

For  the  quadrilateral  element  only,  the  conventional  isoparametric  transfor- 
mation is  used  to  determine  the  strains  of  the  general  element  from  the  shape 
functions  and  derivative  formulas  of  the  parent  element  referenced  to  its 
cartesian  coordinate  system.  This  transformation  does  not  appear  in  Figure  6 
because  it  occurs  at  an  earlier  stage  of  the  calculations. 

Stiffness  Matrix  Transformations 

The  stiffness  matrices  of  the  elements  are  initially  derived  in  the  baseplane 
coordinate  system  and  include  all  freedoms  of  the  element.  For  the  quadri- 
lateral, these  total  58  freedoms,  of  which  40  are  nodal  freedoms  which  are 
retained  for  the  final  problem  solution,  and  18  are  HMN  freedoms,  of  which 
10  are  eliminated  by  explicit  strain  constraints,  4 are  deleted  to  avoid 
probable  inter-element  incompatibilities,  and  4 are  eliminated  by  a minimum 
energy  constraint.  All  constraints  are  on  the  elemental  level.  For  the 
triangular  element,  there  are  a total  of  51  freedoms,  of  which  27  are  nodal 
freedoms  which  are  retained  for  the  final  problem  solution,  and  24  are  HMN 
freedoms,  of  which  6 are  eliminated  by  explicit  strain  constraints,  6 are 
deleted  to  avoid  verified  inter-element  incompatibilities,  and  12  are  elimi- 
nated by  a minimum  energy  constraint.  All  constraints  are  at  the  elemental 
level.  The  stiffness  matrix  transformation  sequences  are  shown  on  Figure  7 
for  the  two  elements.  The  designates  on  the  figure  have  the  following 
meanings: 

This  is  the  transformation  which  imposes  explicit  constraints  on 
the  higher  polynomial  strain  terms.  It  operates  on  the  stiffness 
matrix  as  a conventional  generalized  coordinate  transformation  of 
the  form  C^KC. 

This  is  a simple  removal  of  rows  and  columns  to  eliminate  freedoms 
which  might  cause  inter-element  incompatibil ities. 


HMN 


■ » • 

DELETE 
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MPE  - This  is  a conventional  stationary  potential  energy  reduction 

performed  on  the  elemental  level. 

TSTAR  - This  is  a generalized  coordinate  transformation  performed  to 

change  from  the  BCIZ  "deformational " freedoms  to  freedoms  which 
have  a consistent  sign  convention  and  permit  rigid  body  motions 
of  elements. 

TD  - See  Figure  6.  This  puts  freedoms  into  the  solution  coordinate 

systems. 

MERGE  - This  is  a conventional  merge  of  elemental  freedoms  to  obtain 
structural  equilibrium  equations. 

The  stiffness  matrix  is  formed  at  the  outset  of  each  load  step,  as  indicated 
on  Figure  6.  The  solution  proceeds  with  iterative  corrections  within  the 
step  until  either  convergence  is  obtained  or  an  input  iteration  limit  is 
reached.  In  the  latter  case,  the  stiffness  matrix  is  reformed,  and  all  trans- 
formations performed  again  as  indicated  on  Figure  7.  The  stiffness  matrix 
formation  and  transformation  accounts  for  about  75%  of  the  computational  time 
on  the  small  problems  studied  to  date. 

Loads  Transformations 

The  elements  in  their  present  form  accept  only  nodal  load  inputs.  These  are 
input  in  the  global  system,  and  can  be  used  directly  in  problem  solutions 
after  a coordinate  transformation  using  TCAP  (see  Figure  6).  However,  the 
computation  of  the  residual  loads  for  iteration  is  done  with  all  of  the 
element  freedoms,  and  the  resulting  load  vector  must  be  transformed  and 
merged  through  the  same  steps  used  for  the  stiffness  matrix.  Thus,  Figure  7 
applies  also  to  the  transformation  of  the  elemental  residual  loads.  It  should 
be  noted  that  the  elemental  residual  loads  are  computed  with  reference  to  the 
most  current  updated  element  baseplane.  Hence,  the  TD  matrix  used  (Figure  7) 
must  be  the  most  recent  update  of  that  matrix.  The  MPE  transformation  of  the 
loads  is  the  conventional  one  which  creates  a partial  load  vector  which  is 

saved,  to  be  used  later  in  the  back-substitution  solution  of  the  displacements. 

I Lt  . 
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The  loads  transformations  are  done  for  every  iteration  within  the  load  step,  in 
contrast  with  the  stiffness  matrix  calculations,  which  are  dene  much  lo.'S 
frequently. 


I 
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Displacement  Transformations 

The  incremental  solution  of  a load  step  or  an  iteration  is  initially  referred 
to  the  solution  coordinate  systems  (Figure  5,  6).  A number  of  transformations 
and  other  calculations  are  made  on  these  incremental  values,  as  shown  by  Figure 
8.  The  figure  shows  the  incremental  solution,  consisting  of  nodal  displacement 

quantities  and  called  Aclsoi uti on  coords’  as  the  startin9  Point  of  the  data 
reduction  procedure.  The  calculations  initially  proceed  along  three  separate 
paths:  (1)  TCAP  is  used  to  transform  the  increment  to  global  components  and  to 
update  the  global  coordinates,  element  baseplanes,  and  transformation  matrices 
associated  with  the  baseplanes;  (2)  the  incremental  solution  is  transformed  to 
the  start-of-step  baseplane,  using  the  old  TD  matrix,  and  then  summed  with  the 
prior  accumulated  nodal  displacements  referred  to  this  baseplane;  the  total 
displacements  are  transformed  by  means  of  the  RDOT  matrix  to  obtain  total  nodal 
displacements  referred  to  the  updated  baseplane;  (3)  the  incremental  nodal  dis- 
placements referred  to  the  start-of-step  baseplane  are  used  in  back-substitution 
to  obtain  the  minimum  potential  energy  (MPE)  incremental  freedoms,  Aa,  which  are 
then  summed  to  form  the  running  total  of  these  quantities.  At  this  point  the 
calculation  paths  merge  and  the  explicit  strain  constraints  (HMN  conditions)  are 
imposed,  using  total  rather  than  incremental  strains,  and  directly  computing 
total  HMN  freedom  values.  The  results  at  this  point  include  the  accumulated 
nodal  displacements  q,  referred  to  the  updated  baseplanes,  and  the  total 
accumulated  MPE  and  total  HMN  freedoms,  a and  8.  These  data  suffice  to  compute 
strain,  stress,  and  residual  loads.  The  loads  are  transformed  as  discussed 
above,  merged  to  form  overall  structural  residuals,  and  tested  for  convergence. 
If  convergence  has  not  been  obtained,  a new  increment  of  nodal  displacements, 
designated  as  A(A<l)S0iuti0n  coords’  comPutecl»  using  the  same  solution  coordi- 
nate system  and  stiffness  matrix  as  were  determined  at  the  start  of  the  step. 

If  convergence  has  been  obtained,  the  solution  coordinate  systems  are  updated, 
and  new  TCAP  and  new  stiffness  matrices  are  generated.  In  addition,  a number 
of  data  arrays  are  saved,  as  indicated  on  the  figure.  In  some  cases,  if  con- 
vergence is  slow,  the  stiffness  matrices  are  updated  without  obtaining 
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convergence.  In  this  case,  since  the  residual  loads  remain  referenced  to  the 
old  solution  coordinate  system,  the  latter,  and  TCAP,  are  not  updated.  This 
particular  option  amounts  to  a forced  "yes"  answer  to  the  convergence  test, 
except  that  the  TCAP  is  not  updated.  The  advantage  of  updating  the  stiffness 
matrix  is  due  to  properly  accounting  for  the  effect  of  element  deformed  shape 
on  element  stiffness,  and  to  accounting  for  element  orientation  (updated  base- 
plane)  relative  to  the  solution  coordinate  systems. 

The  HMN  transformation,  as  used  in  stiffness  matrix  transformation,  is  derived 
for  use  by  incrementation.  However,  the  actual  calculation  of  the  HMN  free- 
doms is  a total  rather  than  an  incremental  calculation.  Prior  to  the  data 
processing  of  Figure  8,  HMN  constraint  matrices  in  terms  of  total  rather  than 
incremental  strains  are  formed,  by  simple  changes  in  the  earlier-generated 
incremental  matrices.  Total  HMN  freedoms  are  then  calculated,  avoiding  cumu- 
lative error  in  these  freedoms.  For  the  quadrilateral  only,  the  HMN  matrices 
are  also  updated  prior  to  this  calculation,  to  account  for  the  total  bending 
deformation  accumulated  to  the  current  point  of  the  iterative  calculations. 

The  need  for  this  updating  is  due  to  the  fact  that  the  quadrilateral  basic 
displacement  shapes  are  only  second  degree  forms.  Hence,  the  basic  membrane 
strain  states  are  linear,  and  cannot  compensate  for  nonlinear  effects  due  to 
the  simplest,  i.e.,  constant  curvature  and  constant  twist,  lateral  displace- 
ments. The  triangular  element,  on  the  other  hand,  having  basic  displacement 
forms  which  are  cubic,  can  at  least  partially  compensate  for  nonlinear  effects 
due  to  the  simplest  bending  deformations  with  its  basic  membrane  strain  states. 
The  updating  of  the  HMN  matrices  was  not  necessary  for  the  triangular  element. 

2.3  Solution  Procedures 

Solution  procedure  development  has  required  a large  amount  of  effort  in  the 
present  research,  even  though  the  primary  aim  of  the  work  has  been  element 
evaluation  and  improvement.  The  cause  of  this  is  in  the  nature  of  the  resid- 
ual loads.  In  elements  of  the  types  under  study,  in  which  nonlinear  strains 
due  to  relatively  large  bending  displacements  are  included  during  each  solu- 
tion step,  very  large  membrane  stresses  and  membrane  residual  loads  occur. 

In  combination  with  the  bending  slopes,  these  residuals  create  residual 


bending  loads  and  moments  which  usually  completely  dominate  the  bending 
behavior  of  the  structure.  This  behavior  becomes  more  severe  as  the  plate 
or  shell  becomes  thinner,  corresponding  to  the  increasing  ratio  of  the  mem- 
brane stiffness  to  the  bending  stiffness.  In  general,  the  bending  residual 
loads  are  in  the  proportion  (A/r)  to  the  applied  bending  loads,  where  A is 
the  deflection  magnitude  and  r is  the  section  (or  plate/shell)  radius  of 
gyration.  The  bending  residual  loads  are  distributed  over  the  structure, 
among  the  different  elements,  according  to  the  local  magnitudes  of  the  slopes. 
Hence,  in  performing  a residual  load  iteration,  the  bending  displacement 
adjustment  tends  to  be  much  larger  than  the  initial  displacement,  and  is 
distributed  quite  differently  over  the  structure.  If  the  iterative  adjust- 
ment is  accepted  at  its  computed  magnitude  and  distribution,  generally  a 
grossly  distorted  deformation  state  and  greatly  increased  residual  load  values 
result.  This  type  of  behavior  almost  always  causes  divergence  of  the  solution 
procedure.  This  situation  is  in  marked  contrast  to  the  behavior  which  occurs 
in  analysis  approaches  in  which  residual  loads  are  either  not  computed,  or 
are  computed  by  approximate  methods  which  effectively  omit  the  effects  of  the 
nonlinear  straining  which  occurs  during  the  increment.  In  these  approaches, 
the  residual  loads  are  small,  and  convergence  is  generally  rapidly  achieved.  How- 
ever, for  strongly  nonlinear  problems,  such  approaches  generally  yield  solu- 
tions which  diverge  increasingly  from  the  correct  solution,  as  the  magnitudes 
of  the  nonl inearly-induced  strains  increase. 

The  simplest  way  of  improving  solution  procedure  behavior,  in  fully  nonlinear 
analysis,  is  to  scale  the  magnitude  of  the  iterative  correction  increment. 
Procedures  often  are  coded  with  factors  such  as  0.5,  for  example,  to  be 
applied  to  these  increments.  Though  this  may  at  times  be  successful,  it  is 
an  unsatisfactory  method  in  two  ways:  the  proper  scale  factor  generally  varies 
widely  during  a complete  problem  solution;  the  method  fails  to  address  the 
fact  that  the  iterative  increment  is  distributed  incorrectly  over  the  structure. 
These  two  difficulties  are  the  principal  obstacles  to  achieving  a rapidly  con- 
verging stepwise/iterative  procedure  for  nonlinear  analysis.  It  is  required 
to  achieve  "intelligent,"  programmed  procedures  for  controlling  the  size  of  the 
iterative  increment,  and  also  for  controlling  its  "shape."  The  concept  of 
"shape,"  as  used  here,  includes  both  the  distribution  of  deformations  over  the 
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entire  structure,  and  also  differences  in  deformation  magnitudes  between 
different  types  of  deformation,  e.g.,  membrane  as  opposed  to  bending  motions, 
etc.  The  present  research  has  developed  primitive,  but  effective,  procedures 
for  dealing  with  these  two  difficulties,  and  has  thereby  achieved  convergence 
in  problems  which,  in  earlier,  conventional  versions  of  the  solution  procedure, 
were  not  sol vable. 

Subsequent  paragraphs  will  describe  the  developed  solution  procedure.  However, 
it  is  worthwhile  to  digress  at  this  point  in  order  to  consider  the  implications 
of  this  work  with  regard  to  stepwise-linear  solution  procedures  in  general.  It 
has  been  demonstrated  by  the  solution  procedure  work  of  this  research  that  both 
magnitude  and  direction  are  incorrect  in  stepwise-linear  solutions  of  nonlinear 
problems.  The  errors  begin  with  the  initial  application  of  load,  and  continue 
through  all  iterations  and  subsequent  load  steps.  It  has  also  been  verified 
that  updating  the  stiffness  matrix  tends  to  alleviate  these  difficulties. 
Clearly,  the  basic  source  of  the  error  is  in  the  stiffness  matrix  itself,  both 
as  regards  overall  stiffness  magnitude,  and  as  regards  stiffness  coupling 
effects  between  different  deformation  types,  particularly  between  membrane  and 
bending  displacements.  This  total  problem  would  be  solved  by  using  a nonlinear 
stiffness  description  of  the  structure,  i.e.,  by  using  a stepwise-nonlinear 
approach  for  the  solution  of  strongly  nonlinear  problems.  Such  an  approach  is 
available,  and  is  described  in  References  8 and  9.  It  is  termed  the  "static 
perturbation"  method,  and  achieves  fully  stepwise-nonlinear  performance  without 
a great  deal  more  computational  effort  than  the  conventional  stepwise-linear 
approaches.  For  many  problems,  the  stepwise-nonlinear  approach  would  very 
likely  be  less  expensive  as  well  as  more  reliable,  since  fewer  iterations  and 
fewer  stiffness  matrix  updates  would  generally  be  required. 

Solution  Procedure-Convergence  Acceleration 

The  overall  solution  process  consists  of  a set  of  load  steps,  within  each  of 
which  is  a set  of  iterations,  each  of  which  is  based  on  the  residual  loads 
corresponding  to  the  immediately  previous  displacement  state.  This  latter 
state  may  result  from  either  an  input  load  step  application  or  an  iteration  for 
a residual  load  application.  The  first  significant  improvement  in  convergence 
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resulted  from  performing  iteration  steps  which  alternately  permitted  all  struc- 
tural freedoms  to  respond,  or,  alternatively,  permitted  only  the  membrane  free- 
doms to  respond.  In  this  case,  the  membrane  freedoms  are  those  approximately 
parallel  to  the  plate  or  shell  surface,  as  defined  by  the  solution  coordinate 
systems  (Figure  5).  The  effect  of  the  membrane  freedom  iteration  was  to  allow 
the  structure  to  relieve  much  of  the  nonlinear  strain  and  stress  induced  by 
the  previous  bending  displacement  increment.  This  nonlinear  strain  and  stress 
relief  serves  to  greatly  reduce  the  magnitudes  of  the  residual  loads,  with  the 
result  that  the  following  (all  freedom)  increment  will  have  much  smaller  bend- 
ing displacements,  and,  hence,  much  less  further  nonlinear  strain  and  stress 
generation.  The  purely  membrane  increment  can  be  thought  of  in  two  ways:  as 
a post-increment  correction  of  the  "shape"  error  of  the  previous  all -freedom 
increment;  and  as  a plausible  physical  action,  which  a plate  or  shell  structure 
would  naturally  undertake  to  relieve  equilibrium  imbalances  and  achieve  a 
reduced  potential  energy  level.  These  are  distinct,  but,  of  course,  basically 
equivalent  interpretations.  This  solution  procedure  feature  influenced 
deformation  shape  in  a primary  way,  and,  through  residual  load  reduction, 
influenced  increment  magnitude  secondarily.  Convergence  was  thereby  obtained 
for  problems  which  had  previously  diverged  except  for  very  small  load  increments. 

It  was  found,  however,  that  for  some  problems  the  residual  loads  still  tended 
to  be  large  enough  that  convergence  was  not  obtained  except  for  small  load 
steps.  These  problems  were  generally  those  with  more  elements  or  problems  in 
which  the  overall  structural  stiffness  depends  strongly  on  displacement  magni- 
tude (e . g . , the  un-prestressed  membrane  problem).  The  basic  difficulty  was 
felt  to  be  in  the  magnitude  of  the  displacement  increments,  and  how  they 
were  distributed  over  the  structure. 

To  improve  this  aspect  of  the  convergence,  a procedure  was  implemented  in 
which,  in  the  all-freedom  iterations,  the  amplitude  of  the  increment  was 
arbitrarily  varied  over  a certain  range,  say  50%  to  150%  of  the  computed  value, 
and  a state  of  minimum  residual  was  sought.  For  this  purpose  a measure  of  the 
residual  based  on  a root-sum-square  over  all  the  residuals  was  used.  The 
minimum  was  sought  on  the  basis  of  a quadratic  fit  of  the  residual-measure 
versus  the  amplitude  factor.  This  procedure  performed  well  several  times,  but 
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in  general  had  the  characteristic  of  extremely  slow  convergence.  The  reason 
was  found  to  be  as  follows:  the  residual-measure  for  all  amplitudes  of  the 
all-freedom  increment  (say,  50%,  100%,  and  150%)  generally  exceeded  the  value 
of  the  residual -measure  of  the  previous  membrane-only  increment,  in  rough 
proportion  to  the  magnitude  factor.  Thus,  a minimum  residual  was  not  found. 

The  solution  procedure  was  coded  to  reset  the  range  of  the  search,  to,  say, 

25%,  50%,  75%,  and  repeat  the  calculation.  Generally,  the  failure  to  find  a 
minimum  repeated  itself.  Ultimately,  the  procedure  accepted  a very  low  ampli- 
tude factor  for  the  all-freedom  increment,  and  thus  failed  to  make  appreciable 
progress  toward  convergence. 

The  basic  cause  of  these  difficulties  was  the  "shape"  of  the  all-freedom  incre- 
ment, specifically  the  incorrect  amplitudes  of  its  membrane  as  compared  to  its 
bending  freedoms.  To  remedy  this,  a repeated  use  of  the  membrane-only  incre- 
ment was  implemented  as  follows: 

For  each  amplitude  of  the  all-freedom  increment,  (say,  50%,  100%,  150%), 
residual  loads  are  computed,  and,  based  on  these  new  residuals,  a membrane- 
only  increment  is  computed  and  added  to  the  factored  all-freedom  increment. 

Total  displacements,  residual  loads,  and  the  residual-measure  are  recom- 
puted for  this  "double  increment",  and  this  residual-measure  is  identified 
as  belonging  to  the  particular  amplitude  factor  used. 

The  search  for  a minimum  residual -measure  is  done  using  these  "hybrid" 
residual  load  states  and  residual-measures. 

The  chosen  amplitude  factor  is  used  for  a final,  fourth- time  calculation, 
of  the  displacement  increment  and  the  residual  loads.  This  calculation 
also  uses  a membrane-only  substep. 

In  using  this  procedure,  the  single,  membrane-only  increment  is  no  longer 
necessary  and  it  was  removed  from  the  solution  procedure.  The  double-search 
calculations  are  done  for  each  load  increment  application  as  well  as  for  the 
residual  load  iterations. 
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This  procedure  has  been  successful  on  all  problems  to  which  it  has  been 
applied.  It  is  by  no  means  a fully  optimized  procedure,  and  it  is  clearly 
primitive  and  costly.  Nevertheless,  it  has  provided  convergence  for  previ- 
ously badly  divergent  cases.  It  is  felt  that  the  principal  importance  of 
the  procedure  is  that  it  verifies  that  increment  "shape"  control  and  amplitude 
control  are  both  required,  on  a per-increment  basis,  to  obtain  convergence  to 
a general  class  of  strongly  nonlinear  problems.  This  has  led  to  the  strong 
conviction  that  some  type  of  stepwise  nonlinearity,  based  on  a stiffness 
matrix  which  varies  within  each  increment,  is  essential  to  cost-effective 
solution  procedures  for  strongly  nonlinear  finite  element  analysis. 

Several  other  items  were  found  to  be  necessary  for  convergence,  in  addition 
to  the  overall  method  described  above: 

The  stiffness  matrix  requires  frequent  updating,  particularly  while 
relatively  large  bending  displacement  increments  are  occurring.  This 
greatly  improves  the  quality  of  the  computed  increments. 

The  stiffness  matrix  updating  includes  a re-calculation  of  the  geometric 
stiffness  matrix.  In  this  calculation,  the  stresses  of  the  previous 
converged  steps  are  used  until  the  membrane  stress  state  has  recovered 
from  its  initial  large,  nonlinearity-induced,  excursion.  The  residual- 
measure  magnitude  is  used  to  control  this  decision  process. 

The  residual  measure  has  been  modified  to  consider  only  bending  direction 

(w,  0 , 9 ) residuals.  This  was  found  essential  to  prevent  the  search 
x y 

(optimization)  procedure  from  optimizing  toward  zero  membrane  residuals 
at  the  expense  of  the  bending  residuals.  This  did  occur,  and  caused 
convergence  toward  very  small  displacement  amplitudes,  in  some  cases. 

The  iterative  increment  is  checked  numerically  to  assure  that  it  does  not 
exceed  the  basic  increment  magnitude  of  the  load  step.  If  it  does,  it  is 
scaled  overall,  prior  to  the  search  calculations. 
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If  the  search  calculations  do  not  find  a minimum  residual-measure,  the 
appropriate  end  of  the  search  range  is  accepted  as  a satisfactory  ampli- 
tude factor.  This  has  proven  satisfactory,  and  saves  computation  time 
compared  to  a complete  search  range  change  and  re-calculation. 

Increment  rotation  magnitudes  are  checked,  and  the  increment  is  prevented 
from  rotations  large  enough  to  violate  the  incremental  small  angle 
assumption  (about  20°).  j 

Solution  Procedure--User  Input 

Because  different  types  of  nonlinear  problems  are  best  handled  with  somewhat 
different  calculation  sequences  and  details,  certain  data  can  be  input  by  the 
user,  as  controls  over  the  solution  process. 

The  degree  of  refinement  required  for  converged  solution  states  is  an  input 
item.  Too  small  a tolerance  on  this  value  can  increase  computer  costs  excessively. 

The  number  of  iterations  to  be  performed  before  the  stiffness  matrix  is 
updated  is  controlled  by  user  input.  Thus,  for  example,  a set  of  values  such 
as  1,2, 2, 3 specifies  that,  after  the  load  step  application,  the  stiffness 
matrix  is  updated  immediately.  Thereafter,  two  iterations  are  performed,  the 
stiffness  is  updated,  two  more  iterations  are  performed,  etc. 

The  total  number  of  permitted  stiffness  matrix  updates  is  specified  by  the 
user.  Since  the  updates  are  major  computer  cost  items,  this  provides  user 
control  over  costs  of  solutions  which,  for  some  reason,  are  converging  too 
slowly. 

The  search  range  for  the  "double-search"  procedure  is  controlled  by  user 
input  values  of  the  center  point  and  the  range  on  either  side  over  which 
the  search  for  minimum  residuals  will  be  performed.  Thus,  for  example, 
the  values  of  .80  and  .35,  respectively,  will  cause  the  search  for  the 
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minimum  residual  state  to  take  place  with  amplitude  factors  .45,  .80,  and 
1.15.  This  control  is  particularly  useful  for  structures  which  are  known 
to  stiffen  or  soften  markedly  as  they  accumulate  deformation. 

In  order  to  handle  problems  with  either  "following"  types  of  boundary  con- 
ditions, or  with  conventional  boundary  conditions,  the  user  can  input  par- 
ticular fixed  axes  with  respect  to  which  boundary  conditions  will  be 
imposed.  These  are  referred  to  the  global  system.  Alternatively,  the 
solution  coordinate  system  triads,  which  are  convected  with  the  deforma- 
tion can  be  used  to  enforce  boundary  conditions. 

In  addition  to  the  above,  of  course,  the  user  may  elect  to  insert  zero  load 
steps.  These  steps  cause  a general  cleanup  of  all  coordinate  system  trans- 
formations and  data  updating.  This  has  not  been  found  necessary  in  the  prob- 
lems solved  to  date,  with  the  most  recent  solution  procedure. 

Finally,  the  program  has  been  structured  to  allow  operation  in  the  HMN-mode, 
or  as  a conventional  nonlinear  program  (using  the  basic  AZI  and  BCIZ  elements). 
This  has  permitted  comparisons  of  HMN  and  non-HMN  element  performance. 

2.4  A Consistent  Transformation  for  Finite  Rotational  Freedoms 

With  isoparametric  shell  elements  of  the  AZI  type  (Reference  6),  or  with  corre- 
sponding types  of  curved  beam  elements,  a general  finite  element  model  will 
have  three  translational  and  three  rotational  freedoms  at  each  mode.  The 
three,  rather  than  two,  rotational  freedoms  at  each  node  are  required  in  order 
to  handle  elements  which  intersect  at  other  than  180-degree  angles.  (Three 
rotational  freedoms  are  also  required  if  beam  torsional  behavior  is  to  be 
included. ) 

The  translational  freedoms  may  be  transformed  between  the  solution  and  element- 
baseplane  coordinate  systems,  using  a simple  3x3  matrix  of  direction  cosines, 
and  that  transformation  may  be  applied  exactly  to  finite  incremental  or  cumula- 
tive translations.  In  the  case  of  finite  rotations,  however,  the  resulting 
configuration  is  dependent  upon  the  order  in  which  the  rotations  are  performed. 


This  order  dependence  cannot  be  considered  when  incremental  rotations  are 
calculated  during  the  usual  matrix  solution  procedure  (there  the  rotations 
are  considered  to  be  infinitesimal).  After  a series  of  such  rotation  incre- 
ments are  calculated,  transformed  from  solution  to  baseplane  systems,  and 
summed  to  cumulative  values,  an  inconsistency  develops  relative  to  the 
deformed  configuration  of  adjacent  elements.  That  is,  a cumulative  error  is 
created  which  is  equivalent  to  the  existence  of  physical  gaps  (slope  discon- 
tinuities) between  elements.  This  inconsistency  affects  the  residual-force 
equilibrium  computations,  so  that  it  cannot  be  eliminated  by  the  iteration 
procedure. 

In  order  to  develop  a consistent  transformation  for  the  rotational  freedoms, 
imagine  that  at  each  node  the  following  entities  exist. 

1)  a small,  rigid  globule  of  material 

AAA 

2)  an  arbitrarily  defined  solution-coordinate-system  triad  X-Y-Z 

A 

3)  a unit  normal  vector  NT  for  each  (Ith)  element  joined  to  this  node  (the 

1 A 

initial  orientation  of  Nj  in  the  undeformed  configuration  is  normal  to 
the  Ith  element  baseplane) 

These  three  entities  are  rigidly  attached  together,  so  that  they  translate 
and  rotate  as  a single  rigid  body  during  successive  load  increments.  Figure 
9a  illustrates  the  concept  for  the  simple  case  of  a node  where  three  curved- 
beam  type  elements  are  joined.  The  three  entities  and  the  adjoining  elements 
are  depicted  there  in  the  initial  undeformed  configuration.  Figure  9b  shows 
one  of  the  elements  in  a later  deformed  state,  with  a new  vector  N which  we 

A 

define  as  being  a unit  normal  to  the  deformed  baseplane.  Because  N has  rotated 
with  the  rigid  material  globule  at  the  node,  and  this  rotation  is  in  general 
different  from  that  of  the  element  baseplane,  the  vectors  N and  N are  not 

A 

parallel.  The  difference  between  these  two  vectors  (N-N)  will  provide  a 
convenient  measure  of  the  element  nodal  rotations. 


1; 
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Specifically,  the  consistent  element  rotations  are  determined  according  to 
the  following  steps. 

AAA 

1)  Given  start-of-increment  solution  coordinate  triad  X°-Y°-Z°  at  a node. 

This  triad  can  be  used  as  the  base  vectors  for  expressing  other  vector 
quantities. 

2)  Perform  the  matrix  solution  procedure  to  obtain  nodal  incremental  transla- 

AAA 

tions  and  rotations  (expressed  in  terms  of  X°-Y°-Z0  system  triad). 

3)  Using  the  nodal  translations,  determine  new  element  baseplane  orienta- 
tions and  associated  normal  vectors  N. 

4)  Using  the  nodal  incremental  rotations,  calculate  new  end-of-increment 

AAA 

solution  coordinate  triad  X'-Y'-Z'  at  a node.  (This  calculation  is 
somewhat  arbitrary  due  to  the  order  dependence  of  finite  incremental 
rotations.  However,  the  new  triad  will  be  used  consistently  for  all 
elements  joining  the  node,  so  that  any  error  introduced  will  be  elimi- 
nated by  the  residual  force  iteration. 

A 

5)  Since  N is  rigidly  attached  to  the  solution  triad,  it  is  given  in  terms 
of  X'-Y'-Z',  and  can  then  be  transformed  to  the  base  vectors  X°-Y°-Z°. 

6)  Cumulative  element  nodal  rotations  are  then  computed  from  the  difference 

A 

(N-N),  by  taking  its  vector  dot  product  with  the  element  baseplane 
coordinate  axes. 

This  procedure  allows  the  calculation  of  cumulative  element  nodal  rotations, 
which  have  no  cumulative  error.  That  is,  they  are  consistent  among  all 
elements  adjoining  the  node,  so  that  no  slope  discontinuity  effects  are 
allowed  to  accumulate. 
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2.5  Discussion  of  Numerical  Results 


l I 


The  triangular  and  quadrilateral  HMN  elements  have  been  used  to  solve  a large 
number  of  simple  problems.  The  early  numerical  work  dealt  principally  with 
the  triangular  element,  and  demonstrated  that  the  performance  of  this  element 
was  satisfactory,  both  in  overall  element  linear  and  nonlinear  behavior,  and 
also  in  regard  to  the  HMN-related  behavior.  Example  #5  of  this  section 
illustrates  this  work.  It  was  decided  on  the  basis  of  this  work  and  early 
experience  with  the  quadrilateral  element  that  the  latter  is  the  better 
element  of  the  two  by  a wide  margin.  Consequently,  the  majority  of  HMN- 
element  evaluation  was  done  with  the  quadrilateral.  This  work  is  covered  by 
examples  #1  through  #4  of  this  section.  The  basic  difficulty  with  the  HMN- 
BCIZ  element  is  believed  to  be  related  to  a known  deficiency  of  the  basic  BCIZ 
element,  namely,  that  it  fails  to  maintain  inter-element  slope  conformity. 

The  goals  of  the  numerical  evaluations  are  as  follows: 

To  understand  the  basic  behavior  of  the  elements,  considering  particu- 
larly the  influences  of  coordinate  systems,  the  HMN  function  activity, 
and  the  effects  of  initial  curvature. 


To  compare  the  behavior  of  the  HMN-elements  with  non-HMN-elements, 
considering  both  large  and  small  deflections,  and  both  initially  flat 
and  initially  curved  elements. 


To  develop  solution  procedure  concepts  and  methods  suitable  to  provide 
good  convergence  properties  for  problems  using  HMN-elements. 


These  goals  have  been  met  satisfactorily,  and  a thorough  evaluation  of  the 
elements  has  been  made.  It  has  not  been  possible,  however,  to  demonstrate 
the  elements  in  nonlinear  problems  of  practical  importance  (e.g.,  shell 
buckling),  due  to  the  solution  procedure  difficulties  which  were  encountered 
and  also  the  fact  that  the  computer  programs  were  designed  for  use  with  very 
few  elements.  It  has  also  not  been  undertaken  to  verify  nonlinear  predictions 
against  available  exact  solutions,  despite  the  original  intention  to  make 
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such  checks.  Little  additional  effort  would  be  required  at  this  point,  due  to 
the  convergent  solution  procedures  developed,  to  make  such  checks  for  simple 
problems. 

One  of  the  conclusions  reached  in  the  evaluations  of  the  HMN-AZI  element  is 
the  desirability  of  a formulation  change  to  a modified  type  of  isoparametric 
element  having  displacements  entirely  nodally  defined,  with  explicit  HMN 
constraints  on  higher  degree  strain  polynomials  replaced  by  potential  energy 
constraints.  The  bases  of  this  conclusion  are  contained  in  the  discussions 
of  the  numerical  examples  of  this  section,  particularly  examples  #2  and  #3. 

Example  1 : Pinched  Cylinder 

Figure  10  shows  a cylindrical  shell  pinched  by  line  loads  180°  apart.  To 
study  this  problem,  a quarter-circle  model  using  four  AZI-HMN  elements  was 
analyzed,  using  the  boundary  conditions  indicated.  Both  zero  and  non-zero 
values  of  Poisson's  ratio  were  used.  This  problem  was  a valuable  example  in 
many  respects.  At  the  outset,  it  proved  a severe  test  of  solution  procedures, 
and  led  to  several  improvements  in  these  methods.  In  addition,  the  pinched 
cylinder  problem  has  several  particularly  interesting  features  for  the  present 
investigation:  (1)  it  is  a sensitive  indicator  of  the  differences  between 
shallow  shell  theory  and  deep  shell  theory;  (2)  it  shows  very  strong  effects 
of  the  action  of  the  HMN  freedoms;  (3)  it  illustrates  clearly  the  relationship 
between  the  use  of  the  baseplane  coordinate  system  and  the  HMN  functions  and 
constraints;  (4)  it  illustrates  a difficulty  inherent  in  nonlinear  plate 
bending  analysis  with  nonzero  Poisson's  ratio;  (5)  it  affords  an  opportunity 
to  consider  a case  of  nonlinearity  of  purely  geometric  (negligible  nonlinear 
strain  or  stress  effects)  origin. 

The  radius  of  the  cylinder  is  1"  and  its  width  is  0.4".  Two  thicknesses  were 
used:  .025"  and  .100".  Loads  from  .34#  to  67.9#  were  applied  vertically  at 
the  upper  edge  of  the  quarter-circle  structural  model,  and  both  ends  of  the 
quarter-circle  were  restrained  against  rotation.  Deflections  at  the  load  up 
to  .098"  were  computed.  This  deflection  value  is  9.8%  of  the  radius  and  3.92 
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times  the  thickness  of  .025".  At  this  deflection  a noticeable  softening  of 
the  structure,  due  to  geometry  change,  has  occurred,  as  compared  to  deflections 
for  lower  load  levels.  The  apparent  softening  is  estimated  to  be  15%,  based 
on  the  deflection  under  the  load.  The  non-HMN  solutions  yield  deflections 
ranging  from  atout  1/3  to  1/4  of  the  HMN  values,  depending  on  the  displacement 
amplitude.  A check  of  the  computed  displacements  at  very  small  load  levels 
against  the  comparable  analytical  solution  for  ring  bending  shows  that  the 
four  element  array  used  in  this  problem  is  about  5.7%  too  stiff.  This  is  a 
reasonable  error  for  a four  element  array,  using  constant  curvature  elements, 
for  a problem  of  this  type.  Figure  10  shows  the  deflected  shape  of  the 
structure,  plotted  to  scale,  and  the  force-deflection  curve,  for  the  case  of 
t = .025". 

The  stress  resultants  in  this  problem  are  in  general  quite  rapidly  varying 
over  the  90°  arc,  and  are  strongly  dependent  on  the  number  of  elements  used, 
the  value  of  Poisson's  ratio,  the  shell  thickness  (R/t),  and  the  type  of 
shell  theory  (deep  or  shallow)  used.  This  type  of  complexity  was  found  in 
the  present  study,  and  also  in  Reference  7 for  the  linear  case.  In  general, 
finite  element  internal  stresses  cannot  be  compared  directly  with  applied 
loadings,  except  for  very  small  elements.  Even  in  this  case,  often  the  stress 
resultants  are  difficult  to  interpret,  except  on  the  basis  of  some  type  of 
average  value  over  entire  elements.  The  reason  for  this  lies  in  the  fact  that 
the  finite  element  method  achieves  equilibrium  with  respect  to  generalized 
loads,  which  are  work-equivalents  of  the  stresses,  rather  than  with  respect  to 
actual  boundary  values  of  the  stresses  themselves.  In  the  present  problem, 
evaluation  of  the  stress  resultants  was  therefore  made  qualitatively.  Stress 
resultant  values  were  compared  for  the  HMN  and  non-HMN  cases  and  for  deep 
versus  shallow  shell  strain  equations  (see  Section  2.1,  Shell  Strain  Displace- 
ment Equations).  It  was  found  that  the  non-HMN  stress  predictions  in  all 
cases  showed  large  oscillation  behavior  both  within  and  between  elements, 
while  the  HMN  stress  results  were  very  smooth  and  plausible  in  character  over 
the  entire  structure.  The  membrane  stress  resultant  parallel  to  the  circular 
arc  was  found  to  be  sensitive  in  both  magnitude  and  distribution  to  whether 
the  deep  or  shallow  type  of  shell  equations  were  used,  for  the  case  of  the 
thick  (R/t  = 10)  shell.  The  other  stress  resultants  were  not  influenced  by  the 
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choice  of  strain  equations.  Most  of  the  stress  resultants  were  strongly  and 
plausibly  affected  by  Poisson's  ratio;  this  is  discussed  later  in  this  section. 

The  extremely  strong  influence  of  the  HMN  freedoms  and  constraints  can  be 
understood  quite  clearly  in  this  problem.  It  is  recalled  that  the  element 
displacement  states,  referred  to  the  updated  element  baseplane,  are  cartesian. 
Considering  an  initially  curved  element,  such  as  those  employed  in  the  present 
example,  it  is  easily  seen  that  a flattening  of  the  element,  due  to  bending 
displacements  only,  will  cause  coireciable  compressive  straining  near  its 
ends.  The  compressive  strain  due  to  sjch  a flattening  is  given  by 
w°’x*w,x’  w^ere  w°  *he  ''''’itial  curve  of  the  element,  w is  the  incremental 
elastic  deformation,  and  x is  the  baseplane  cartesian  coordinate  along  the 
length  of  the  element.  The  HMN  displacements  were  specifically  constructed 
to  eliminate  this  type  of  straining,  which  has  a second  degree  behavior  in  x, 
and,  through  the  HMN  constraint  equations,  this  strain  and  the  corresponding 
stress  are  eliminated.  The  important  thing  to  note  is  that  this  effect,  as 
described  here,  is  purely  linear.  The  high  degree  strain  is  due  to  linear 
behavior  in  the  presence  of  initial  element  curvature,  and  results  from  the 
use  of  the  cartesian  baseplane  coordinate  system.  The  continued  deformation 
of  the  element  of  course  causes  the  development  of  nonlinear  strains  of  the 
same  type,  which  are  also  removed  from  the  deformation  state  b the  HMN 
functions  and  constraints.  It  might  be  argued  that  a change  fv’om  the  cartesian 
coordinates  to  the  shell  mid-surface  curvilinear  coordinates  would  eliminate 
the  high  polynomial  degree  strains  in  the  linear  case,  this  removing  the  need 
for  the  HMN  functions  and  constraints.  While  this  is  true,  of  course,  it 
would  in  turn  cause  much  more  serious  calculation  problems  in  the  nonlinear 
case.  This  is  discussed  under  "Shell  Equations,"  and  "Coordinate  Systems  and 
Updating,"  in  Section  2.1. 

Several  plate  bending  problems,  including  the  present  one,  have  showed  a 
nonlinear  effect  involving  bending  in  the  presence  of  a non-zero  Poisson's 
ratio.  Consider  a plate  or  shell  element  initially  flat  in  the  y-direction, 
but  either  curved  or  flat  in  the  x-di recti on.  When  the  element  is  bent  into 
an  elastic  curve  in  the  x-direction,  because  of  Poisson's  ratio  it  attempts 
to  respond  by  becoming  concave  in  the  y-direction  on  its  surface  of  bending 


tensile  stress.  This  action  causes  slope  changes  of  the  plate  or  shell, 
given  by  w,^,  which  create  y-direction  nonlinear  membrane  tensile  strains 
proportional  to  (w,y)2,  particularly  near  the  edges  of  the  element.  The  HMN 
elements  eliminate  these  strains  through  the  HMN  v-direction  membrane  displace- 
ment functions.  The  non-HMN  elements  cannot  modify  these  strains,  and  the 
corresponding  stresses,  in  a similar  way,  because  the  basic  AZI-element  membrane 
functions  do  not  provide  displacements  of  suitable  polynomial  degree.  Hence, 
the  non-HMN  elements  experience  a strong,  nonl inearly-induced  resistance 
against  the  cross-curvature  due  to  Poisson's  ratio,  and  respond  by  severely 
limiting  the  amplitude  of  the  y-direction  curve  of  the  elements.  In  the 
present  problem,  the  HMN  element  under  the  load  experiences  a lateral  curve 
of  amplitude  .0012",  while  the  non-HMN  element  amplitude  is  .00028,  for  the 
case  of  the  67.9#  load  and  the  .025"  thickness.  Due  to  the  lateral  curve, 
the  x-direction  membrane  stress  for  the  HMN  element  is  significantly  affected 
due  to  the  resulting  geometry  change.  Near  the  load,  the  edges  of  the  elements 
tend  to  move  toward  the  center  of  the  arc  as  a result  of  the  Poisson's  ratio- 
induced  lateral  curvature.  This  causes  a rather  large  x-direction  compression 
stress  near  the  element  edges,  and  a large  tension  stress  near  the  central 
portion  of  the  element.  This  behavior  is  undetectable  for  the  non-HMN  case, 
for  two  reasons:  the  lateral  curvature  is  negligible;  the  x-direction  membrane 
stress  is  very  badly  behaved  because  of  the  lack  of  the  HMN  function  action 
to  handle  the  initial  curvature-induced  stresses  (see  discussion  above). 

Reference  7 shows  a strong  dependence  of  numerical  results  for  this  problem, 
for  R/t  = 10  and  v = 0,  between  finite  element  solutions  based  on  deep  shell 
and  shallow  shell  equations.  This  effect  was  evaluated  with  the  present 
element,  also  for  R/t  = 10  and  v = 0.  The  results  indicated  that  the  deflec- 
tions and  all  stress  resultants  except  the  circumferential  membrane  stress 
are  only  slightly  affected  by  the  inclusion  or  omission  in  the  strain- 
displacement  equations  of  the  terms  such  as  w°,x*u,x  and  w»x*u>x»  ln 
circumferential  bending  curvature.  Moreover,  the  bending  displacements  in 
both  cases  check  closely  with  the  exact  solution  (linear).  In  the  present 
formulation,  based  on  cartesian  displacement  definitions,  including  or  omitting 
these  types  of  terms  is  the  only  way  in  which  shallow  or  deep  shell  types  of 
theories  can  be  simulated.  It  is  concluded  that,  for  the  type  of  strain 


displacement  equations  used  in  this  research,  the  usual  errors  associated 
with  shallow  theory  do  not  occur  to  any  noticeable  extent.  This  appears  to 
be  due  to  the  fact  that,  for  the  cartesian  displacements  used,  the  usual 
approximations  of  the  shallow  shell  type  no  lonqer  cause  any  significant 
approximation  in  the  equations  for  the  curvatures  and  the  twist.  It  was 
noted  that  the  circumferential  stress  resultant  was  quite  strongly  affected 
by  the  change  of  equations.  This  may  or  may  not  be  significant,  however, 
because  this  particular  resultant  is  not  well  predicted  for  the  coarse 
finite  element  model  used,  and  is  extremely  sensitive  to  the  shell  equations 
used  in  a way  which  is  dependent  on  element  size  (Reference  7). 

Example  2:  Cantilever  Plate  with  Initial  Curvature 

Figure  11  shows  the  problem  under  consideration.  This  is  the  interesting 
"carpenters  tape"  problem,  in  which  nonlinearity  occurs  because  of  flattening 
of  the  cross-section.  Many  related  problems  are  of  importance  in  structural 
analysis.  This  structure  acts  initially  as  a simple  cantilever  beam.  For 
either  up-loads  or  down-loads,  the  edges  of  the  beam  tend  to  deflect  in  such 
a way  as  to  flatten  the  initial  lateral  curve  of  the  cross-section,  thereby 
reducing  its  effective  bending  moment  of  inertia  and  causing  loss  of  stiffness. 
The  reasons  for  solving  the  problem  in  the  present  work  were:  to  determine 
whether  the  HMN  functions  would  act  (as  they  should)  to  facilitate  the  flat- 
tening of  the  section,  which  for  the  non-HMN  case  should  be  greatly  attenuated; 
to  study  the  HMN-element  behavior  for  a case  of  initial  curvature  transverse 
to  the  principal  loading  direction.  This  problem  illustrated  several  very 
interesting  facets  of  the  behavior  of  HMN  elements. 

Overall,  the  solutions  obtained  showed  noticeable  loss  of  stiffness,  due  to 
cross-sectional  flattening,  for  the  HMN-element  cases,  and  a slight  tendency 
for  stiffening  with  the  non-HMN  element.  The  cross-section  showed  the  expected 
flattening  for  the  HMN-element  cases  and  also  for  the  non-HMN  element  for  the 
up-load  case  only.  For  the  down  load  cases,  the  non-HMN  element  showed  a 
reverse-flattening,  i.e.,  an  increased  "cupping"  of  the  section.  The  nonlinear 
effects  were  small  at  the  load  levels  studied,  and  the  HMN  and  non-HMN  results 
were  comparable  in  overall  stiffness.  Flattening  of  the  section  reached 


37 


about  10«  of  the  initial  curved  shape,  and  end  deflections  reached  about 
.36",  which  is  about  22%  of  the  span.  This  deflection  is  much  greater  than 
either  the  thickness  or  the  initial  out-of-flatness  of  the  section. 

The  flattening  behavior  computed  with  the  HMN  elements  was  well  behaved  and 
completely  plausible.  The  somewhat  greater  stiffness  and  reduced  flattening 
for  the  up-load  case  is  consistent  with  the  easily  observed  behavior  of  a 
carpenter's  tape  subjected  to  upward  and  downward  directions  of  loading. 

This  behavior  is  also  consistent  with  the  general  effect  of  cross-curvature 
induced  by  Poisson's  ratio,  and  may  possibly  be  related  to  this  effect.  The 
ability  of  the  HMN-elements  to  flatten  is  due  directly  to  the  action  of  the 
HMN  functions,  since  it  is  through  these  functions  that  the  development  of 
large,  y-direction  membrane  strains,  due  to  the  flattening,  is  avoided. 

The  anomolous  flattening  action  of  the  non-HMN  elements  is  caused  by  a very 
interesting  and  unexpected  behavior.  Due  to  the  overall  bending  loading  in 
the  case  of  an  upward  load,  the  cross-section  experiences  tensile  stress  at 
the  edges  and  compressive  stress  at  the  crown.  The  reverse  occurs  for  a 
downward  load.  Considering  first  the  downward  load,  it  is  seen  that  Poisson's 
ratio  will  tend  to  create  compressive  y-direction  stresses  near  the  edges  of 
the  element  and  tensile  stresses  near  the  crown.  This  y-direction  stress  is 
roughly  quadratic  in  the  y-coordinate.  At  this  point  a significant  charac- 
teristic of  the  AZI  element,  and  of  many  multi -mode  finite  elements,  is  noted: 
the  element  has  no  facility  at  all,  within  linear  theory,  to  rid  itself  of 
these  particular  Poisson's  ratio-induced  stresses,  because  of  the  character  of 
its  available  displacement  functions.  In  the  present  case,  however,  where 
nonlinear  strains  are  included,  the  element  can  can  rid  itself  of  these 
stresses  almost  completely.  The  means  of  doing  this  is  by  "cupping" 
the  section  to  an  increased  lateral  curvature.  This  action  produces,  through 
the  terms  w°»y*w’y  1n  the  y-direction  membrane  strain,  a tensioning  of  the 
outer  edge  zone,  and,  in  conjunction  with  the  basic  element  v displace- 
ments, a compression  in  the  crown  region.  This  action  tends  to  reduce 
the  potential  energy.  Thus,  for  the  case  of  downward  loading,  the  section 
becomes  more  laterally  curved,  and  the  overall  structural  behavior  becomes  somewhat 
stiffer.  Of  course,  the  well  known  tendency  of  the  section  to  flatten  (because 
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of  equilibrium  effects)  is  present  in  the  non-HMN  as  well  as  the  HMN-element. 

It  simply  occurs  in  the  present  case  that  the  above-described  behavior  due  to 
the  Poisson's  ratio  effect  dominates  the  behavior.  The  flattening  tendency 
is  driven  by  relatively  small  forces,  and  cannot  compete  with  the  very  dominant 
effects  associated  with  the  membrane  energy. 

For  the  upward  load  case,  the  non-HMN  element  must  flatten  its  section  to 
conform  to  the  Poisson's-ratio-related  behavior  as  discussed  above.  In  this 
case  the  Poisson's-ratio-induced  behavior  reinforces  the  natural  tendency  of 
the  section  to  flatten.  The  data  (Figure  11)  show,  however,  that,  while  the 
section  does  flatten  in  this  case,  it  does  so  only  to  about  the  same  degree 
that  the  "cupping"  occurred  for  the  downward  load  case.  This  behavior  appears 
surprising,  but  is  easily  explained.  The  explanation  is  again  in  the  dominant 
effect  of  the  membrane  energy.  The  lateral  bending  of  the  section  has  been 
seen  to  create  large  y-direction  membrane  strains,  and  thus  affects  the 
membrane  energy  directly.  For  the  non-HMN  case,  which  cannot  relieve  these 
membrane  stresses  and  strains,  the  equilibrium-driven  flattening  tendency  (at 
small  loads)  is  simply  not  strong  enough  to  overcome  the  inherent,  membrane 
strain-induced  stiffness  against  lateral  bending  (flattening).  Thus,  the 
non-HMN  element  solution  shows  a flattening  or  cupping  of  the  section  which 
is  totally  dominated  by,  and  essentially  serves  to  exactly  balance,  the 
Poisson's  ratio-induced  y-direction  stresses.  As  observed,  the  effect  is 
essentially  exactly  reversed  for  the  upload  and  download  cases.  For  the  HMN- 
element,  however,  the  lateral  bending  is  only  resisted  by  the  small  plate 
bending  stiffness,  because  the  HMN  functions  eliminate  any  y-direction  membrane 
stress  participation.  Thus,  the  HMN  elements  show  section  flattening  in 
direct  response  to  the  normal  flattening  tendency,  as  driven  by  the  equations 
of  equilibrium  in  the  deformed  structural  shape. 

One  final  observation  concerning  the  HMN-element  results  for  this  problem  is 
noted.  The  Poisson's  ratio-induced  y-direction  membrane  stresses  due  to  the 
initial  bending  stress  distrubution,  seen  to  dominate  the  non-HMN  element 
behavior,  gives  rise  to  absolutely  no  response  in  the  HMN-element.  The 
element  does  not  have  any  first-order  ability  within  its  basic  functions  to 
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rid  itself  of  these  stresses,  because  of  the  simple  forms  of  the  available  v 
displacement  functions  (nor  does  the  non-HMN  element).  It  cannot  use  the  HMN 
functions  either,  because  the  HMN  constraints,  which  are  the  only  means  of 
activating  the  HMN  functions,  respond  to  strain  only,  and  not  to  stress. 
Finally,  it  cannot  rid  itself  of  these  stresses  through  lateral  curvature,  as 
the  non-HMN  elements  do,  because  the  HMN  constraints  prevent  the  development 
of  y-direction  deformations  in  this  case  (of  course,  this  particular  behavior 
is  very  undesirable  anyway).  The  conclusion  which  results  from  this  discussion 
is  that  the  HMN  element  conceptual  basis  should  be  extended  to  apply  to  high 
polynomial  degree  stresses  as  well  as  strains.  In  effect,  this  would  require 
committing  the  entire  HMN  process  to  the  control  of  the  potential  energy 
theorem. 

Example  3:  Cantilever  Plate  with  End  Restraint 

Figure  12  shows  the  problem  under  consideration.  This  structure  is  highly 
nonlinear  due  to  the  effect  of  the  displacement  constraint  at  the  right  end. 

Its  deflection  magnitudes  are  controlled  almost  entirely  by  the  large  tension 
stresses  induced  by  this  constraint,  which  are  distributed  almost  uniformly 
over  the  span.  The  HMN  and  non-HMN  elements  produce  virtually  identical 
results  under  these  circumstances  (about  1%  difference).  The  problem  was 
solved  primarily  as  a test  of  the  solution  procedure.  In  order  to  obtain 
convergence  in  this  case,  it  is  necessary  to  avoid,  or  otherwise  attenuate 
the  effect  of,  the  extremely  large  residual  loads  which  result  from  the  initial 
steps  of  the  solution,  during  which  the  stiffening  effect  due  to  the  end 


condition  is  not  present  in  the  stiffness  matrix.  The  present  solution 
procedure  (Section  2.3)  handles  this  by  looking  at  a range  of  amplitudes  of 
the  initial  solution  step  displacements.  A state  of  minimum  residual  is 
sought  within  this  range  to  determine  the  optimal  amplitude.  In  the  present 
case,  a minimum  does  not  occur  within  the  search  range,  and  its  lower  limit 
amplitude  is  selected.  The  second  step  proceeds  similarly.  Because  of  the 
rapid  change  of  stiffness  with  deflection,  the  stiffness  matrix  is  updated 
frequently.  This  type  of  solution  procedure  avoids  the  requirement  that  the 
user  set  intelligently  chosen  amplitude  factors  for  the  attenuation  of  the 
residual  load  iterative  procedure.  It  also  avoids  divergence  even  with 


relatively  large  load  steps.  Prior  to  implementation  of  the  optimum-amplitude 

| search-type  procedure,  solution  attempts  diverged  for  this  problem. 

j The  figure  shows  a plot  of  the  deflected  shape  of  the  edge  of  the  plate.  The 

centerline  deflections  are  slightly  different  due  to  the  .3  value  for  Poisson's 
ratio.  The  initial  linear  steps  had  an  end  deflection  of  .155"  for  zero 
Poisson's  ratio  and  .148  for  v = .3.  The  linear  theory  end  deflection  for 
zero  Poisson's  ratio  is  .157".  The  straight-line  character  of  the  deflection 
shape  for  the  outer  two  elements  is  a result  of  the  nonlinearity-induced 
tensile  stresses.  The  nonlinear  solution  gives  about  20 % of  the  linear 
solution  value  of  the  end  deflection,  at  maximum  load. 

Example  4:  Twisting  of  a Square  Plate  Fixed  on  One  Side 

Figure  13  illustrates  the  problem,  which  was  solved  with  both  a single  element 
and  with  four  elements,  as  indicated.  For  the  four  element  model  the  applied 
loads  are  distributed  over  the  end  such  that  the  deflections  there  are  linear 
in  y.  The  applied  torque  is  250  inch-pounds.  The  character  of  this  problem 
is  dominated  by  the  use  of  the  total  fixity  condition  at  the  supported  edge. 

The  problem  is  significantly  nonlinear,  with  the  initial  step,  linear  solution 
maximum  deflections  exceeding  the  final  converged  deflections  in  the  ratios 
of  1.65  for  the  single  element  and  1.87  for  the  four  element  case.  The 
membrane  stresses  are  (and  should  be)  quite  large  for  this  problem,  and  are 
important  influences  on  the  torsion-bending  behavior.  Because  of  the  fixed 
edge  and  the  fact  that  the  structure  would  "prefer"  to  adopt  a deformation 
state  of  constant,  pure  twist,  rather  than  the  torsion-bending  state  demanded 
by  the  edge  constraint,  the  single  element  model  is  much  too  stiff  in  the 
linear  case.  The  linear  first  step  single  element  maximum  displacement  is 
.161",  while  for  the  four  element  model  the  value  is  .224".  The  deflection 
data  are  tabulated  on  the  figure. 

For  the  case  of  the  fixed  edge,  this  problem  is  one  in  which  the  HMN  elements 
as  formulated  in  this  research  demonstrate  unacceptable  behavior.  The  diffi- 
culty is  in  the  x-direction  variation  of  the  y-direction  direct  strain,  ey. 

This  particular  behavior  (cy  vs.  x)  is  not  considered  by  the  HMN  constraints 
or  the  HMN  supplementary  displacement  functions.  The  slope  w,y  is  a second 


degree  function  of  x because  of  the  fixed  edge.  Thus,  the  nonlinear  behavior 
generates  a fourth  degree  function  of  x in  the  strain  ey.  No  displacement 
options  are  available  to  counteract  this  strain,  which  remains  in  the  deforma- 
tion state  as  a source  of  excessive  energy  and,  hence,  excessive  stiffness. 
Solutions  of  this  problem  with  both  the  single  and  four  element  models  show 
that  the  stress  varies  as  a high  degree  (4th)  function  of  x.  Since  this 
occurs  for  both  the  v = 0 and  v = .3  cases,  it  is  clearly  caused  by  slope  w,y  and 
the  inability  of  the  present  HMN  formulation  to  deal  with  the  problem  in  a 
successful  manner.  The  largest  values  occur  at  the  loaded  end  of  the 
plate,  of  course,  because  the  lateral  slope  is  largest  there. 

The  figure  shows  a "carpet"  plot  of  the  variation  of  vs  y on  five  x=constant 
lines  on  the  plate.  Both  the  single  element  and  four  element  cases  are 
shown.  For  the  single  element  model,  for  which  v = 0,  oy  is  constant  in  y,  and 
the  implied  cross-plots  would  show  to  be  quartic  in  x.  For  the  four 

element  model,  for  which  v = .3,  the  effect  of  Poisson's  ratio  is  to  cause 
to  be  quadratic  in  y over  each  individual  element.  The  distribution  is 
symmetric  about  the  x-centerline,  but  a small  discontinuity  in  a occurs 
between  the  elements  which  join  each  other  at  x = .5.  Again,  the  implied 
cross-plot  variation  of  with  x is  quart ic.  The  four  element  model  shows 

comparable  values  to  the  single  element  model  overall,  but  has  considerably 
larger  values  at  the  loaded  end.  The  cause  of  these  larger  stresses  lies 
partly  in  the  added  twist  angle  of  the  four  element  model  (about  23%),  and 
partly  in  the  rather  sharply  peaked  distribution  of  with  y.  The  latter  is 
probably  due  to  the  effect  of  Poisson's  ratio. 

The  figure  indicates  the  regions  of  tension  and  compression  a . In  both 
models  this  stress  is  compressive  over  the  central  region  (0<X<1;  .25<y<.75) 
and  tensile  over  the  edge  regions  (0<X<1;  o<y<.25;  .75<y<l).  The  cause  of 
this  behavior  is  of  course  the  tendency  of  the  edge  zones  to  foreshorten,  due 
to  their  slope,  which  tendency  is  resisted  by  the  central  zone  of  the 
plate.  This  gives  rise  to  a rather  large  level  of  membrane  shear  stress 
as  well.  The  single  element  and  four  element  models  have  roughly  comparable 
maximum  values  of  the  x-direction  direct  stress  resultant  and  the  shear 
stress  resultant,  as  follows:  four  element  model,  Nx  = 21300,  Nx  = -8000, 
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N = 5800;  single  element  model,  N = 7300,  N = -5000,  N = 8300.  The 
xy  x x xy 

larger  values  of  N^  for  the  finer  model  reflect  its  ability  to  adequately 
represent  the  x-direction  distribution  of  the  stress;  the  Nx  versus  x plot 
shows  a maximum  at  x = .5. 

The  HMN  functions  have  the  potential  to  cause  inter-element  incompatibility 
of  the  x-direction  displacement  between  the  joining  elements  at  x = .5.  The 
data  were  checked  for  the  finer  model,  and  this  incompatibility  was  found  to 
be  present,  though  small.  This  results  from  the  strong  twisting  action  of 
this  problem,  resulting  in  large,  nonlinear  membrane  shear  stresses.  This  is 
considered  to  be  an  unsatisfactory  behavior  on  the  part  of  the  HMN  elements 
in  their  present  formulation. 

The  conclusions  reached  from  this  analysis  of  nonlinear  torsion-bending  are: 

the  HMN  formulation  must  deal  with  the  nonlinear  strains  c , e , and  e as 

x y xy 

a coupled  set,  no  one  of  which  can  be  optimized  in  any  way  independent  of  the 
others;  the  possibility  of  inter-element  incompatibility  should  be  eliminated 
by  controlling  the  HMN  displacements  by  nodal  freedoms.  These  conclusions 
point  toward  a formulation  in  which  the  HMN  theory  is  based  on  potential 
energy  constraints  only,  and  manipulates  displacements  which  are  entirely 
nodally  controlled. 

Example  5:  BC1Z-HMN  Cantilever  Beam 

Figure  14  shows  a three  element  cantilever  beam  bent  by  end  loads.  The 
initial  problem  solutions  with  this  model  used  z-direction  concentrated 
loadings  at  nodes  4 and  5.  The  results  for  this  loading  showed  the  beam  to 
be  about  20%  too  stiff,  and  also  yielded  a free  end  moment  of  45%  of  the 
fixed  ended  moment,  within  the  outer  pair  of  elements.  The  net  end  moment  at 
the  free  end  was  zero,  of  course,  as  the  result  of  the  combined  action  of  all 
three  elements.  The  residual  x-direction  stresses  along  the  outer  edges  of 
the  element  were  proportional  to  the  square  of  the  load  (hence  related  to 
nonlinear  behavior),  about  1%  of  the  potential  nonlinear  stress,  and  of 
quadratic  variation  along  the  edge,  rather  than  the  expected  linear  variation. 

These  results  indicate  that  though  the  HMN  procedure  is  functioning  correctly,  some 
unexpected  nonlinear  effect  is  preventing  the  basic  BCIZ  membrane  functions 
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from  eliminating  the  nonlinear  x-direction  stresses.  It  was  concluded  that, 
probably  due  to  the  BCIZ  element's  characteristic  of  permitting  inter-element 
slope  incompatibilities,  the  minimum  energy  state  of  the  element  actually 
retains  small  membrane  stress  levels.  More  important,  the  20%  stiffness 
error  and  the  erroneous  bending  moment  distribution  were  judged  unacceptable, 
even  for  the  coarse  idealization  employed,  because  these  elements  use  cubic 
displacement  functions,  and  should  have  been  very  accurate  for  this  problem. 
Therefore,  this  problem  was  investigated  further  under  different  loadings,  as 
indicated  in  Figure  14.  The  loadings  are,  respectively,  concentrated  nodal 
moments,  generalized  loadings  based  on  a distributed  edge  moment  along  edge 
4-5,  and  loadings  deduced  directly  from  the  element  stiffness  matrices  to 
obtain  a constant  curvature  state.  The  numerical  results  obtained,  as  given 
in  the  figures,  show  the  model  to  be  accurate  within  1%  for  the  nodal  moment 
case,  to  be  about  20%  soft  for  the  generalized  (work-equivalent)  loading,  and 
to  be  about  6%  soft  for  the  loading  derived  from  the  stiffness  matrices. 

These  percentages  are  based  on  "proportional"  loading,  and  refer  to  the 
magnitudes  of  the  moment  loadings  themselves.  The  accuracy  obtained  for  the 
pure  nodal  moment  case  is  excellent,  but  of  questionable  meaning.  That 
obtained  in  the  generalized  load  case  is  clearly  unacceptable.  The  6%  error 
for  the  constant  curvature  case  is  especially  surprising,  since  the  elements 
should  have  responded  to  this  deformation  in  an  exact  way.  It  is  noted  that 
the  residual  stresses  are  identically  zero  for  this  case.  In  the  case  of  a 
structure  incorporating  many  of  these  elements,  it  is  difficult  to  say  what 
sort  of  accuracy  would  be  obtained.  For  constant  curvature,  the  BCIZ  elements 
should  not  exhibit  inter-element  slope  incompatibilities.  It  is  concluded 
that  further  investigation  along  these  lines  would  be  required,  before  the 
value  of  the  BCIZ  element  as  a basis  for  a stability  element  could  be  estab- 
lished. Because  of  the  strong  coupling  between  the  HMN  freedoms  and  the 
slopes,  the  lack  of  slope  continuity  inherent  with  the  BCIZ  element  may  be  a 
serious  difficulty.  These  difficulties  caused  the  BCIZ-HMN  element  to  be 
dropped  from  the  present  research,  in  favor  of  the  AZI-HMN  element. 
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2.6  Alternatives  to  the  HMN- Element  Approach 

This  section  attempts  to  discuss  alternative  methods  which  may  provide  equiva- 
lent gains  to  those  obtained  with  the  HMN  approach,  in  application  to  strongly 
nonlinear  problems.  The  discussion  derives  from  basic  reasoning  concerning  the 
goals  and  difficulties  of  solving  nonlinear  problems,  and  not  from  direct 
comparative  numerical  evaluations,  except  for  the  specific  AZI-element  based 
HMN  versus  non-HMN  comparisons  of  the  previous  section.  Thus,  no  claim  of 
completeness  is  made  here.  In  addition,  the  matter  of  HMN-element  performance 
of  initially  curved  elements  in  linear  analysis,  which  was  discussed  in  Sections 
2.1  and  2.5,  is  not  specifically  addressed  again  here.  However,  most  of  the 
conclusions  drawn  herein  for  nonlinear  analysis  are  also  applicable  to  the 
linear  analysis  of  initially  curved  elements,  for  element  types  whose  formula- 
tions are  specifically  tailored  for  nonlinear  problems. 

The  accurate  and  cost-effective  finite  element  solution  of  nonlinear  problems 
would  appear  to  depend  on  such  factors  as  the  following: 

1)  ease  of  discretization  and  reasonable  control  of  problem  size  (degrees  of 
freedom)--the  use  of  reasonably  large  and  reasonably  uniform  sizes  of 
elements 

2)  avoidance  of  excessive  computer  run  times  and  costs,  as  regards  element- 
related  calculations  (stiffness  matrices,  decompositions,  residual  loads) 

3)  reliable  and  reasonably  fast  convergence  of  residual  load  iterations--the 
ability  to  use  relatively  large  load  steps 

4)  avoidance  of  the  likelihood  of  an  aborted  solution  or  a solution  which  has 
converged  to  an  incorrect  result 

The  HMN  procedure  was  aimed  principally  at  item  1);  i.e.,  it  is  designed  to 
solve  linear  and  nonlinear  problems  with  equal  accuracy  for  the  same  element 
sizes.  At  this  point  it  cannot  be  concluded  whether  it  is  particularly  advan- 
tageous as  regards  items  2)  and  3).  It  does  well  with  respect  to  item  4), 
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converging  to  correct  results  in  cases  where  non-HMN  elements  of  the  same  size 
and  basic  formulation  do  not. 

It  appears  that  all  of  the  technical  gains  of  the  HMN  elements  can  be  obtained 
with  conventional  elements  by  simply  using  smaller  element  sizes  in  a fully 
nonlinear  solution.  This  is  a very  attractive  alternative  (although  admittedly 
most  currently  available  conventional  elements  do  not  retain  full  nonlinearity 
of  strains),  and  one  which  probably  would  not  increase  computer  costs  exces- 
sively for  most  problems.  For  very  large  problems,  however,  the  rapidly  increas- 
ing cost  of  decomposition  of  the  structural  stiffness  matrix,  as  the  number  of 
elements  is  increased,  would  become  a serious  problem.  Perhaps  the  most  serious 
objection  to  this  approach  is  that  the  analyst's  task  is  made  considerably  more 
difficult  in  this  case.  The  necessary  element  sizes  are  not  easy  to  guess 
before  the  first  numerical  solutions  are  obtained,  and  it  will  be  true  in  many 
problems  that  markedly  different  sizes  of  elements  will  be  required  in  different 
areas  of  a structure,  due  solely  to  the  local  degree  of  nonlinearity  of  behavior. 
For  the  HMN  elements,  on  the  other  hand,  the  analyst  is  free  to  use  the  same 
discretization  which  would  be  suitable  for  linear  problems.  In  addition,  once 
a solution  has  been  obtained  with  small,  fully  nonlinear,  conventional  elements, 
the  question  of  whether  a significant  nonlinearity-induced  error  has  occurred 
would  have  to  be  addressed.  A direct  approach  to  answering  this  question  would 
be  to  compare  the  membrane  strain  magnitudes  caused  by  the  linear  and  the 
nonlinear  deformations.  This  is  not  always  a satisfactory  approach,  however. 

The  difficulty  is  that  the  conventional  elements  are  constrained  (by  the  global 
equations  and  the  residual  loads)  such  that  they  underestimate  the  bending 
curvatures,  particularly  in  zones  of  significant  nonlinearity.  Thus,  the  true 
magnitudes  of  the  nonlinear  strains  may  be  considerably  larger  than  the  apparent 
values  given  by  the  solution  obtained.  Based  on  these  considerations,  it  would 
appear  that,  overall,  the  HMN  type  of  element  has  a significant  advantage  over 
the  use  of  smaller  conventional  elements,  in  most  problems. 

Another  possible  approach  involves  limiting  the  types  of  terms  which  are  included 
in  the  calculation  of  nonlinear  strains.  For  example,  in  the  case  of  an  element 
which  retains  linearly  varying  membrane  displacements,  the  computation  of 
nonlinear  strains  might  consider  only  the  average  bending  slopes  of  the  element 
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in  the  calculation  of  the  direct  strains.  A comparable  assumption  would  be 
made  for  the  membrane  shear  strain.  By  this  means  the  nonlinear  strains  are 
kept  to  the  same  functional  forms  as  the  linear  strains,  and  hence  excessively 
stiff  behavior  is  not  caused  specifically  by  higher  degree  strain  polynomials 
due  to  the  nonl inearity.  This  particular  approach  has  a technical  drawback, 
however:  the  elements  will  not  respond  to  the  tendency  of  tensile  stresses  to 
straighten  an  element,  or  of  curvature  increases  to  shorten  an  element  (or  the 
reverse  for  compression,  curvature  decrease).  This  drawback  limits  these 
elements  to  use  in  small  element  sizes  for  nonlinear  problems.  Since  the  type 
of  error  involved  here  acts  directly  on  the  high  energy  overall  tension/ 
compression/shear  membrane  stress  state  of  the  element,  rather  than  on  only  a 
higher  degree  polynomial  component  of  these  stresses,  it  would  appear  to  have  the 
potential  for  causing  serious  errors  in  equilibrium.  It  is  noted  that  taking 
this  type  of  liberty  with  the  strain  equations  necessarily  affects  the  coupling 
between  nonlinear  membrane  direct  strain  and  shear  strain.  If  this 
approach  were  used  for  an  element  such  as  the  AZI  element,  in  which  the  membrane 
displacements  have  a second  degree  polynomial  capability,  and  hence  a first 
degree  membrane  stress  capability  (second  degree  for  shear),  there  is  not  avail- 
able an  efficient  truncation  of  the  bending  slope  polynomial.  Use  of  the 
average  slope  only  would  appear  to  be  an  overly  severe  approximation,  but 
retention  of  the  linear  term  in  the  slope  already  exceeds  the  basic  element 
membrane  strain  polynomial  capability.  Thus,  this  basic  approach  to  handling 
nonlinear  behavior  would  appear  to  be  most  efficient  in  application  to  the 
simplest  type  of  element  as  regards  membrane  displacement--that  in  which  the 
membrane  displacements  are  linear  functions.  Such  elements,  however,  are 
usually  not  accurate  for  either  linear  or  nonlinear  analysis  of  curved  shell 
elements  except  in  very  small  element  sizes.  This  is  because,  for  curved 
shells,  generally  membrane  displacement  behavior  is  as  rapidly  varying  as 
bending  behavior,  except  in  local  edge-bending  zones.  In  addition,  without  the 
use  of  competent  membrane  functions,  there  is  often  difficulty  in  providing 
finite  elements  with  strain-free  rigid  body  motions,  especially  for  large 
deflection  analysis. 

A similar  approach  to  the  above  would  be  one  in  which  the  fully  nonlinear 
strain  is  computed,  but  the  result  is  filtered  in  some  way,  retaining  only 
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lower  (basic  element)  degree  polynomials  for  stiffness  matrix  and  residual  load 
calculation.  We  are  not  aware  of  such  an  approach  being  used,  and  it  is  diffi- 
cult to  guess  what  sort  of  performance  it  would  provide.  It  would  appear  that 
the  method  has  a basic  flaw  in  that  the  membrane  strain  state  and  the  bending 
slope  state  do  not  properly  correspond  to  each  other.  Thus,  bending  displace- 
ments might  be  incorrectly  predicted,  because  their  membrane  (high  energy) 
consequences  are  not  fully  accounted  for.  This  is  much  the  same  sort. of  objec- 
tion as  the  arguments  given  above  against  the  use  of  average  slopes  only  in  the 
nonlinear  strain  equations. 

It  would  be  possible  to  use  a very  complex  element  as  regards  all  types  of 
freedoms.  Thus,  for  example,  an  isoparametric  element  with  5 nodes  per  side 
might  be  used.  This  would  provide  adequate  membrane  strain  capability  to 
handle  nonlinear  strains  caused  by  bending  slopes  corresponding  to  a constant 
curvature  state.  Of  course,  the  possibility  of  very  high  degree  nonlinear 
strains  exists  in  this  case,  because  of  the  high  degree  bending  displacements 
employed.  The  constant  curvature  terms  are  dominant,  however.  It  would  be 
possible  to  use  relatively  large  elements  of  this  type  in  nonlinear  problem 
solutions,  with  results  comparable  to  those  obtained  with  the  HMN  type  of 
element.  This  appears  to  be  a superior  approach,  and  is  similar  to  the  recom- 
mended approach  arrived  at  on  the  basis  of  the  present  research. 

All  of  the  above  possibilities  relate  to  problem  solutions  in  which  nonlinear 
strains  are  computed,  and  the  resulting  nonlinear  stresses  affect  the  residual 
loads  and  stiffness  matrices  of  the  elements.  There  are  also  methods  which 
perform  stepwise  solutions  in  which  nonlinear  effects  are  totally  ignored,  so 
far  as  explicit  calculation  is  concerned.  These  methods  always  use  a "geometric" 
stiffness  matrix  to  represent  the  effort  of  structural  shape  change  on  the 
equilibrium  equations.  This  is  probably  the  most  common  type  of  method  employed 
for  solving  nonlinear  problems  in  current  finite  element  programs.  The  method 
commonly  takes  two  forms:  in  one  case,  no  iteration  is  done,  and  the  structural 
geometry  is  updated,  forming  elements  of  new  shapes,  at  the  end  of  each  step; 
in  the  other  case  a single  step  is  done  and  iterated,  with  the  iteration  based 
not  on  residual'  loads  but  on  re-forming  the  "geometric"  stiffness  matrix  to 
reflect  end-of-step  values  of  the  stresses.  In  the  first  case,  nonlinearity  is 


48 


handled  implicitly  by  the  updating  of  the  shapes  of  the  elements.  However,  the 
nonlinearity  which  occurs  during  a step  is  not  included  because  residuals  are 
not  evaluated.  This  approach  is  known  to  yield  problem  solutions  which  deviate 
increasingly  and  significantly  from  correct  solutions,  even  if  very  small  steps 
are  used.  In  the  second  method,  no  nonlinearity  is  considered  at  all.  At  the 
completion  of  a single  stepwise  calculation,  the  end-of-step  stress  state  is 
used  to  evaluate  an  updated  stiffness  matrix,  and  the  stepwise  calculation  is 
repeated.  The  end-of-step  stresses  are  computed  using  linear  theory.  This 
method  obtains  a correct  solution  only  if  nonlinear  strain  effects  are  negligibl 
It  would  be  possible  to  combine  the  methods,  obtaining  a procedure  which  should 
deviate  more  slowly  from  the  exact  solution  that  the  above-described  conven- 
tional stepwise  method.  In  developing  such  a method  for  nonlinear  problems,  it 
would  be  necessary  to  deal  carefully  with  the  interpretation  of  the  relationship 
between  element  stresses  and  forces,  due  to  changing  element  shapes.  This 
would  prove  somewhat  complex,  and  it  would  appear  that  recourse  to  a fully 
nonlinear  approach  might  be  nearly  as  easy  to  develop.  Nevertheless,  this  type 
of  quasi-nonl inear  analysis  may  have  merit  in  some  types  of  problems. 
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FIGURE  7:  SEQUENCE  OF  STIFFNESS  MATRIX  TRANSFORMATIONS  FOR  STABILITY  ELEMENTS 
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3.0  SUMMARY  AND  CONCLUSIONS 


This  report  concludes  a study  of  finite  element  analysis  of  plate  and  shell 
problems  with  strong  geometric  nonlinearities.  The  specific  area  of  study  is 
the  influence  of  the  nonlinear  strains,  which  generally  are  of  complex  distribu- 
tions over  the  elements,  on  element  stiffnesses  and  accuracy  of  problem  solu- 
tions. The  means  of  investigation  was  through  newly  developed  special  types  of 
finite  elements,  in  which  the  conventional  displacement  functions  are  supple- 
mented by  added  membrane  displacements  of  higher  order  forms.  The  added  func- 
tions are  used  to  modify  the  nonlinear  membrane  strains  in  order  to  obtain  the 
simple  strain  distributions  characteristic  of  linear  analysis.  A combination 
of  constraints  is  used  to  accomplish  this.  They  include  explicit  constraints 
which  set  to  zero  the  higher  order  polynomial  functions  in  the  nonlinear  strains. 
Conventional  potential  energy-based  constraints  are  also  used.  All  constraints 
are  done  at  the  elemental  level. 

In  order  to  implement  the  elements,  solution  procedures  were  developed,  based 
on  a Lagrangian  formulation,  in  which  element- foil  owing,  updated  coordinate 
systems  are  used  to  allow  total  strain  to  be  computed  without  incrementation. 

The  solution  procedure  uses  linearized  stepping  plus  iteration,  with  added 
features  which  account  for  in-step  nonlinear  behavior. 

Evaluation  of  the  special  nonlinear  elements  was  made  through  comparative 
solutions  obtained  with  these  elements  and  with  conventional  elements  of  the 
same  basic  formulation,  all  using  the  same  solution  procedure.  This  work  has 
led  to  the  conclusions  outlined  in  the  following  paragraphs. 

For  strongly  nonlinear  problems,  and  also  for  linear  problems  with  initially 
curved  elements,  the  special  elements  perform  much  better  than  conventional 
elements  when  relatively  large  element  sizes  are  used.  The  error  characteristic 
of  the  conventional  elements  is  excessive  stiffness,  amounting  to  a factor  of 
two  to  four  in  some  of  the  problems  solved.  The  conventional  and  special 
elements  become  essentially  equivalent  in  performance  for  sufficiently  small 
element  sizes. 


For  nonlinear  analysis  in  which  a fully  nonlinear  strain  representation  is 
included,  the  use  of  linearized  steps  generally  yields  displacement  increments 
of  incorrect  "shape,"  in  the  sense  that  the  relative  membrane  and  bending 
displacements  are  incorrectly  proportioned.  In  addition,  for  many  problems, 
the  error  in  displacement  "shape"  and  the  use  of  stepwise  linearity  causes  the 
displacement  increments  to  be  of  incorrect  amplitude  overall.  These  effects 
result  in  excessive  residual  loads  and  difficulty  in  convergence  of  the  itera- 
tion process.  This  research  has  demonstrated  the  beneficial  effects  of  includ- 
ing in-step  nonlinearity  by  means  of  calculations  made  with  a special  solution 
procedure  which  approximates  the  in-step  nonlinearity.  It  is  recommended  that 
a stepwise  nonlinear  approach  such  as  the  static  perturbation  procedure  be  used 
with  the  special  types  of  nonlinear  elements  under  study  in  this  research. 

The  use  of  a stepping  procedure  which  combines  updating  of  element  baseplane 
coordinate  systems,  total  Lagrangian  strain  calculation,  and  locally  cartesian 
displacement  reference  systems  is  particularly  convenient  and  accurate  for 
nonlinear  analysis  of  shells. 

The  special  nonlinear  element  formulation  used  in  this  research  requires  modifi- 
cation in  order  to  achieve  general  applicability.  The  explicit  constraint  of 
higher  order  strain  polynomials  is  not  a satisfactory  approach,  because  it 
interferes  with  the  proper  coupling  between  the  three  membrane  strains.  It  is 
recommended  that  the  formulation  be  modified,  for  isoparametric  elements,  as 
follows:  (1)  the  supplementary  high  order  membrane  displacement  functions 
should  be  controlled  by  nodal  freedoms,  employing  additional  nodes  for  this 
purpose;  (2)  the  added  nodes  need  not,  and  probably  should  not,  be  used  for 
added  bending  freedoms;  (3)  all  constraints  relative  to  the  supplementary 
functions  should  be  enforced  globally,  based  on  conventional  potential  energy 
principles.  A formulation  of  this  type  would  be  simpler  to  develop  and  imple- 
ment than  the  one  used  in  the  present  research,  and  would  provide  all  of  the 
benefits  demonstrated  for  the  special  nonlinear  elements  herein. 
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